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ABSTRACT 


Some of the largest Martian erosive features were influenced by groundwater, and 
include valley networks, outflow channels, and possibly landslides. 

We argue that hydrothermal systems attending crustal formation processes were 
able to drive sufficient groundwater to the surface to form the Noachian southern 
highlands valley networks, which show a spatial correlation to crustal magnetic 
anomalies, also results of crustal formation. Hydrothermal activity is quantified 
through numerical simulations of convection in a porous medium due to the presence 
of a hot intruded magma chamber. The parameter space includes magma chamber 
depth, volume, aspect ratio, and host rock permeability and porosity. For 
permeabilities as low as 10' 17 m 2 and intrusion volumes as low as 50 km J , the total 
discharge due to intrusions building that part of the southern highlands crust 
associated with magnetic anomalies spans a comparable range as the inferred 
discharge from the overlying valley networks. 

The Hesperian circum-Chryse outflow channels are further manifestations of 
groundwater discharge and Clifford and Parker (2001) suggest that the large volumes 
of water required for their formation flows beneath a confining cryosphere from the 
South Pole where meltwater beneath an ice cap recharges a global aquifer. We argue 
that recharge occurs instead over the nearby Tharsis aquifer at high obliquity, assisted 
by cryosphere melting due to volcanic activity. Numerical simulations quantify the 
strength and duration of outflow discharge given either South Polar or Tharsis 
recharge. The contribution of South Pole recharge given Clifford and Parker (2001) 
aquifer properties is negligible compared to that of the initial Tharsis inventory. 
Tharsis recharge, despite the restrictions of improved aquifer properties, makes a 
significant contribution and, unlike South Pole recharge under the same conditions, 
fulfills discharge requirements. 

Groundwater may have influenced long run-out landslide formation in the Valles 
Marineris. We present simulations of Martian, terrestrial, and lunar landslides that 
gauge the role of pore fluid pressure in reproducing accurate geometries and run-out 
with frictional, Bingham, and fluidization rheologies. The results indicate that pore 
fluid is a necessary component of Martian landslide formation and we suggest 
scenarios that might explain its presence. 


This report is the Ph.D. Thesis of Keith Harrison, who performed the work under 
the supervision of the PI, Dr. Robert Grimm. 
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CHAPTER 1 
INTRODUCTION 

1.1 Observational History of Water on Mars 

Introduction of the specific landforms studied in this work is best begun with a 
brief overview of the observational history of water on Mars. The first discovery of 
its presence was through ground-based, infra-red detection of atmospheric water 
vapor in 1963 (Spinrad et al.). Disappointing prospects for surface water followed 
two years later with dry, lunar-like images photographed from Mars orbit by Mariner 
4, followed by much the same from subsequent Mariner 6 and 7 craft (Leighton et al., 
1965; Leighton et al., 1969; Collins, 1971). Only with images from the first Mars 
orbiter (Mariner 9), and in the late 1 970’s with Viking, were fluvial landforms first 
observed, and more detailed questions concerning surface and groundwater raised. 
The Mars Global Surveyor Mission (1997 to the present), with its global capture of 
high-resolution photographic and topographic data, has not only improved our 
understanding of landforms discovered by Viking, but has brought to light another set 
of smaller, more subtle water-related features, including the geologically young spur- 
and-gully formations (Malin and Edgett, 2000). Finally, over the last two years, the 
Mars 2001 Odyssey Mission has detected (with its Gamma Ray Spectrometer) 
significant quantities of ice in the mid- to high-latitude regolith (Boynton et al., 

2002). Our initial impression of Mars as a dry, lunar-like planet has evolved to one of 
an environment affected over most of its history by the action of water. In particular, 
crustal groundwater has played an important erosive role, and We address questions 
regarding three of its manifestations: valley networks, landslides, and outflow 



channels, which We now introduce individually in the order they appear in later 
chapters. 
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1.2 Valley Networks 

The valley networks are the most common drainage feature on Mars. Distributed 
mostly on ancient Southern Highlands terrain within 65° of the equator and at all 
longitudes, they embody the strongest evidence that a warm climate conducive to 
globally distributed flow of surface water once existed on Mars. Groundwater 
sapping is the predominant process responsible for valley network formation (favored 
over precipitation; Baker, 1992; and references therein) and gives rise to long, narrow 
“U” shaped valleys that often follow zones of weakness, such as fault lines. While 
local patterns of groundwater flow leading to sapping are fairly well understood 
(Dunne, 1980), it is not clear how Martian valleys were supplied with the necessary 
water, a question we address in this work. 

1.3 Long Runout Landslides 

The Valles Marineris make up a 4000 km long system of equatorial grabens that 
appear to have formed in response to tectonic forces in the nearby Tharsis rise 
(Blasius et al., 1977), although recent work suggests the influence of local volcanic 
intrusions (McKenzie and Nimmo, 1998). Failure of the canyon wall rock has 
produced numerous landslide events, the largest of which have left runout deposits on 
the canyon floor up to 80 km long. Long run-outs are not predicted by models of 
simple frictional sliding and on the earth usually involve the action of water pressure 
in the pore spaces of the failed material. The presence of groundwater in the Valles 
Marineris slides was first suggested by Lucchitta (1987), and agrees with evidence of 
water-related processes in the formation of other Valles Marineris landforms, such as 



the interior layered deposits (e.g. McCauley, 1978; Nedell et ai, 1987). In this work, 
we attempt to constrain the requirements on pore fluid pressure in Martian long run- 
out landslide formation. 
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1.4 Outflow channels 

While the valley networks reveal information concerning the global distribution 
of Noachian surface water on Mars, the outflow channels provide insight into the 
volume of the global water inventory later in Martian history, during the Hesperian. 
These large channels, unlike the valley networks, have distinct bedforms indicative of 
water flow (Mars Channel Working Group, 1983) and experienced higher, probably 
catastrophic, flow rates (Baker and Milton, 1974). They are thought to have expelled 
at least 6 x 10 6 km 3 of water into the low-lying northern plains (Carr, 1996). This 
water came from the crustal aquifer, probably reaching the surface through a 
disruption in an overlying icy permafrost zone, but details of the necessarily distal 
sourcing of such large quantities of water remain unanswered. Clifford and Parker 
(2001) put forward a sophisticated model involving sublimation from a northern 
lowlands ocean, deposition and accumulation of a South Polar ice cap, and the 
charging of a global aquifer by water melted at its base, we suggest an alternative, 
regional method of aquifer recharge through nearby high elevation Tharsis during 
periods of high planetary obliquity. 


1.5 Tharsis 

Valley networks, long run-out landslides, and outflow channels are related 
through the influence of the Tharsis rise. Trends in valley network orientation 
contribute to the recent consensus that Tharsis formed during the Noachian (Phillips 
etal ., 2001; Anderson et ciL, 2001). The volcanically active nature of its 
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development suggests a strong effect on climate (Jakosky and Phillips, 2001), 
changing global surface stability and atmospheric content of water and therefore the 
conditions under which valley network erosion took place. The location, morphology, 
and orientation of the Valles Marineris suggest a strong correlation to tectonic forces 
(Blasius et al., 1977) in Tharsis that may have been responsible for triggering 
landslides (Caruso, 2003). Likewise, Tharsis influenced the outflow channels, 
possibly disrupting the cryosphere with thermal perturbations, and providing water in 
necessary quantities with the required hydraulic head. 
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CHAPTER 2 

HYDROTHERMAL SYSTEMS AND VALLEY NETWORKS - NOACHIAN 

GROUNDWATER FLOW 

2.1 Introduction 

Valley networks in the ancient terrain of the Martian Southern Highlands 
represent some of the oldest and most widespread evidence of water-related erosional 
processes on Mars. They occur at most longitudes and at low to mid-latitudes, and 
their time of formation is mostly Noachian, suggesting that the early Martian climate 
might have been warm enough to ensure the stability of surface water. The 
morphology of some valleys indicates that they were fed by groundwater discharge 
(see Baker et al , 1992, Baker, 2001, and Carr, 1996, for reviews; Malin and Carr, 
1999; Malin and Edgett, 2000). A question that follows naturally from this discovery 
is what drove groundwater to the surface, especially on the large spatial scale 
observed. Hydrothermal activity may have played a primary role here. 

Hydrothermal systems are important providers of energy in many terrestrial 
igneous, metamorphic, and sedimentary environments, and have profound 
geochemical and biological implications. They accelerate the cooling of magmatic 
bodies in systems ranging from divergent plate boundaries to individual volcanoes, 
and commonly produce substantial surface discharge of water and steam in the form 
of hot springs, geysers, and submarine vents (black smokers). 

On Mars, evidence remains of past hydrothermal activity (see Farmer (1996), for 
a review). Wilhelms and Baldwin (1989) note the widespread presence of ancient 
volcanic structures in the southern highlands, as suggested by outcrop structures 
along the global crustal dichotomy and in crater floors and other depressions. Fluvial 
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channels are frequently associated with volcanic activity (Brakenridge, 1990; Baker 
et al., 1992) and associated structures that include (but are not limited to) ancient 
volcanoes (Gulick, 1998). Hydrothermal circulation may also have altered the martian 
crust and further produced weathering products and soil (Griffith and Shock, 1997; 
Newsom et al., 1999). The greatest consequence of hydrothermal activity on Mars 
may be its ability to sustain life (Shock, 1996; Farmer, 1996). 

we present here numerical models of the thermal convection of groundwater in a 
porous host rock due to the presence of an intruded magma chamber. An extensive 
portion of the available parameter space is explored in order to quantify the effects 
that magma chamber volume, depth, aspect ratio, host rock permeability and porosity 
have on surface discharge. The results of this general approach are then applied to 
some specific issues. First, we use geochemical constraints to bound the permeability 
of the martian crust. Second, we test the hypothesis that hydrothermal circulation can 
explain the putative correlation observed between valley networks and magnetic 
anomalies (Jakosky and Phillips, 2001). we suggest that large amounts of water were 
circulated throughout the southern-highlands crust due to magmatic intrusion and that 
the portion of this water discharged to the surface can quantitatively account for the 
valley networks preserved since the end of crustal formation. 

2.2 Model 

Two-dimensional, axisymmetric representations of hydrothermal circulation in a 
magma chamber and its host rock are modeled with the USGS code HYDROTHERM 
(Hayba and Ingebritsen, 1994, 1997) and visualized with the help of its graphical user 
interface HTpost (P. S. Hsieh, USGS, unpublished material, 2000). HYDROTHERM 
can model temperatures from about 0 to 1200 °C and pressures from 0.5 to 10,000 
bars, and keeps track of both liquid and gaseous phases of pure water. It solves mass. 
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momentum, and energy balance equations expressed in terms of dependent variables 
pressure and enthalpy. The choice of enthalpy over temperature allows the 
thermodynamic state of the fluid to be specified uniquely under two-phase conditions. 
In some circumstances, the combination of pressure and enthalpy is supercritical. 

The resulting fluid, which technically is neither a liquid nor a gas, is treated by 
HYDROTHERM as an equal mixture of each. In my models, supercritical conditions 
are found almost exclusively within the magma chamber, which at these high 
temperatures is ductile and thus impermeable, obviating the need for a detailed model 
of supercritical flow. Viscosity and density for a particular temperature and pressure 
(in any phase) are obtained from a look-up table. 

The momentum balance equation is Darcy's Law, which for a single fluid phase / 
is 

v,=-— (Vp + ^.gV z) (2.1) 

M, 


where v, is the Darcy velocity, k is intrinsic permeability, p, is dynamic viscosity,/? is 
pressure, p, is density, g is gravitational acceleration, and z is depth. The relative 
permeability k' re , quantifies the reduction of the flow of phase / due to the presence of 
the other phase. The mass-balance (continuity) equation for phase / is: 

l(nS i p i ) + V.(p l y i ) = 0. (2.2) 

ot 

where n is porosity, t is time, and S, is volumetric saturation (S wa ier + S s , eum = I , i.e. the 
medium is fully saturated). The continuity equation for the entire system is the sum 
of the equations for each phase. The energy-balance equation for the entire system is 
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d_ 

dt 


(\-n)p r h r +n^S,p,h, 


+ V 


Z* r ,M V , 


VK m VT = 0 


(2-3) 


where h is enthalpy, T is temperature, K m is medium thermal conductivity, and 
subscript r refers to rock matrix properties. HYDROTHERM solves the equations by 
performing Newton-Raphson iterations on an equivalent finite difference system 
(with the horizontal dimension expressed in radial coordinates) until mass and energy 
residuals fall below specified maximum values. 

The radial extent of the host rock is 20 times that of the magma chamber and the 
vertical extent is 20 km, both of which are sufficient to accommodate flow from all 
magma chambers studied. Fluid is allowed to cross the upper horizontal and right 
vertical boundaries, while temperature and pressure are fixed. The right vertical 
boundary is not intended to represent the limit of a horizontally bounded water 
source, but to provide enough space for a realistic response to a local temperature 
perturbation. As much fluid flows through this boundary as is necessary to balance 
the net flow through the top horizontal boundary. The suitability of the chosen radial 
extent was tested by a baseline model measuring 40 chamber radii across, which 
yielded almost the identical discharge to the original model. Recharge into the model 
is typically less than 10% of discharge, indicating that the strength of discharge is not 
offset significantly by infiltration of runoff. Adaptive boundary conditions 
preventing infiltration for a dry Mars could therefore be neglected. 

The initial pressure distribution is hydrostatic, with a surface value of 1 bar. A 
small geothermal gradient (0.5 °C/km), applied to ensure stable decay of the surface 
discharge, does not otherwise affect the model (the effects of substantial geothermal 
gradients were explored and the results are described below). Numerically stable 
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solutions in this model were possible only with temperatures above 10 °C, and this 
was used for the surface boundary condition. Both the surface temperature and 
pressure were adopted for numerical convenience and do not affect the results, which 
are determined almost exclusively by hydrothermal driving forces, nor are they 
intended to represent Earth-like climatic conditions. The left vertical and lower 
horizontal boundaries are no-flow, and the temperature and pressure are free to vary. 
Grid spacing is approximately 100 m in both directions near the magma chamber and 
at greater horizontal distances increases logarithmically. 

Present martian conditions may include a permafrost layer which must first be 
melted before discharge is produced. Gulick (1998) estimated that the time required 
to melt a 2 km thick layer of permafrost was much less than the lifetime of the 
hydrothermal system, making it unlikely that quantities integrated over the lifetime of 
the system, such as total discharge, should be significantly altered. During the late 
Noachian, when the hydrothermal systems proposed here were active, the permafrost 
layer was likely to be much thinner than 2 km, making melting times even shorter. A 
HYDROTHERM simulation that quantifies the role of ice in our models is described 
below. A more detailed explanation of cryosphere growth in the Hesperian appears in 
the second main section of this work (Chapter 3). 

Our baseline model contains a 50 km J chamber emplaced at a depth of 2 km 
below the surface, into host rock of permeability 10‘ 16 m 2 . The dimensions of the 
chamber are most usefully expressed in terms of its aspect ratio, defined here as the 
ratio of diameter to height ( D/H ), which in the baseline model has a value of 2. The 
volume of the chamber is taken after that modeled by Hayba and Ingebritsen (1997) 
and is of a similar size to magma chambers found under mid-oceanic ridges (Burnett 
et al., 1989). Note that our baseline model is intended only as an example model 
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from which we later deviate extensively, and does not necessarily represent any sort 
of “ideal” martian parameter values. 

The magma chamber is emplaced instantaneously at a temperature of 900 °C. 
Because this approach ignores discharge produced during the finite intrusion process 
and subsequent supersolidus cooling, it produces relatively conservative results. The 
chamber is initially impermeable, but as it cools through a brittle-ductile transition 
(BDT) between 400 °C and 360 °C, it is assumed to fracture and become as 
permeable as the surrounding host rock (Hayba and Ingebritsen, 1997). A semilog 
form is adopted for this transition, wherein the log of the permeability scales linearly 
with temperature. Note that the “permeability” BDT may differ somewhat from the 
classic “deformational” BDT (e.g., Kohlstedt et al., 1995) brought on by overburden 
pressure. All models use a rock density of 2500 kg.m'\ a thermal conductivity of 2.0 
W.nT'.K' 1 , and a porosity of 1%. Deviations from the baseline model include host 
rock permeabilities of 1 O' 17 and 10' 15 m 2 , a porosity of 25%, magma chamber aspect 
ratios of 0.2 and 20, depths of 8.5 and 1 5 km, and volumes of 100 and 2000 km 3 . 

Steam and water fluxes, while fundamental to discharge calculations, may also be 
used to test hypotheses regarding the possible geochemical alteration in the system. 
For each time step, a measure of water-to-rock ratio for reactions above a particular 
temperature threshold may be calculated by measuring the total mass of water that 
passes through the region warmer than the threshold value, and dividing by the 
volume of the region. 


2.3 Results 

In the baseline model (Figure 2.1) an initial peak in the surface discharge 
occurring a few hundred years following magma emplacement is due to thermal 
pressurization (Delaney, 1982). Its peak value is not significantly greater than 
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discharges that occur later in the model, and the short-lived nature of this effect, 
which is less extreme for higher host rock permeabilities, contributes negligibly to the 
total, time-integrated mass of water produced by the system (total discharge), which 
we assume to be of primary importance in valley erosion. Additionally, thermal 
pressurization may be weaker in the more realistic case of a finite duration intrusion 
process. 



Figure 2.1 Surface discharge from the baseline model (magma chamber depth 2 km, 
volume 50 km3, aspect ratio 2, host rock permeability 10‘ 16 m 2 ). The regions 
separated by solid vertical lines in the enlarged portion of the figure denote the 
different phases present within in the model. The magma chamber becomes 
permeable over the interval 0.7 to 57 kyr. The noise on the curve is due mainly to 
discretization effects in the presence of steam, which has much greater velocity than 
the liquid phase. 


The broad peak at 45 kyr (1 ky = 1 ,000 yr) is due to thermal convection of 
groundwater. At first, the chamber is supercritical and impermeable, allowing 
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convection in the surrounding host rock only. Surface discharge peaks when the 
magma chamber has cooled sufficiently to become permeable and admit advection. 
When only the liquid phase remains, discharge becomes smoother, decaying 
exponentially and ceasing at 2.8 Myr. 

2.3.1 Permeability 

For simplicity, permeability was homogeneous in all models. Gulick (1998) used 
a value of 10‘ n m 2 for martian hydrothermal systems in analogy with young, near- 
surface Hawaiian volcanics. However, Ingebritsen and Sanford (1998) report that 
permeabilities in the east rift zone at Kilauea, while high near the surface (10‘ 10 to 10‘ 
9 m 2 ), are much lower at depths of just 1 to 2 km (ltf 16 to 10' 15 m 2 ) in rock of the 
same composition. Manning and Ingebritsen (1999) estimate values of between 10* 17 
and 10’ 14 m 2 for the mean continental permeability between 1 and 10 km depth, and 
we therefore use values from 10‘ 17 to 10' 15 m 2 (Figure 2.2). The upper limit is for 
computational convenience; we demonstrate below that results for higher 
permeabilities can be extrapolated from this range. In simulations involving outflow 
channels (Chapter 3) we exploit the full depth dependence of the Manning and 
Ingebritsen (1999) permeability since the flow in these models responds to more 
subtle topographic driving forces. 

Discharge from the low (10' !7 m 2 ) permeability model has a characteristically 
large thermal pressurization peak followed by a weak main peak. Conduction is the 
dominant mode of heat transfer in this model and the spatially integrated surface 
discharge at any given time (instantaneous discharge) is approximately a tenth that of 
the baseline model. While convection is dominant in the high (10* 15 m 2 ) permeability 
model, convection and conduction play comparative roles in the baseline model. 
Discharge in all three models lasts between 2.5 and 3.0 Myr, but the rate of decay is. 
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in general, proportional to permeability, with the 10' 17 m 2 model falling to a tenth of 
its peak value in 0.20 Myr, the baseline model in 0.18 Myr, and the 10‘ 15 m 2 model in 
0.14 Myr. 



time [kyr] 

Figure 2.2 Surface discharge from a 2 km deep magma chamber of volume 50 krrf 
for host rock permeabilities of 10' 17 , 10' 16 , and 10' 15 m 2 . The total discharge of the k 
- 1CT 16 m 2 model, in kg, appears above its curve. The numbers above the other curves 
denote the fraction of this discharge that their corresponding models produced. 


An important feature of the high permeability model is the presence of a steam- 
dominated zone above the magma chamber during the first several thousand years. 
This phenomenon may have implications for chemical alteration, as described below. 

The relationship between discharge and permeability was investigated using the 
high permeability results of Gulick (1998) as a starting point. She modeled a 10' n m 2 
hydrothermal system which included the magma chamber implicitly through a heat 
flux boundary condition based on the analytical solution to the conductive heat flow 
through the wall of an infinitely long cylinder. Consequently, heat loss through the 
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roof and floor of the chamber was not simulated, and no BDT was encountered on 
cooling. Adopting these limitations in HYDROTHERM by removing the host rock 
above and below the magma chamber affords direct comparison with the results of 
Gulick (1998), including her high permeabilities (Figure 2.3) which are more easily 
simulated with this simplified geometry. Agreement is good, although small 
differences result due to our more realistic magma chamber heat flux, and possibly 
other factors such as numerical solution. Overall, we find that the relationship 
between total discharge and permeability can be represented approximately by 


D - a - be~ cK 


(2.4) 


where D and K are the base 1 0 logarithms of total discharge and permeability 
respectively, and a, b, and c are positive constants (with best-fit values of 38, 0.089, 
and 0.12 respectively). As permeability increases, discharge approaches 
asymptotically a maximum value of a, a trend explained by the fixed initial supply of 
heat in the magma chamber of which a growing, yet ultimately limited fraction 
energizes groundwater convection. For the range of permeabilities covered with our 
own, more complete geometry, D is approximately linear in K. Note that the total 
discharges (whose relative magnitudes are given above each non-baseline model 
curve in Figure 2.2) reflect this linear behavior. 
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Figure 2.3 Discharge results from HYDROTHERM models of a similar geometry to 
those of Gulick (1998). An analytic fit to these data, with its large permeability 
asymptote, is shown. The discharge obtained by the single permeability modeled by 
Gulick is provided for comparison, as are the total discharges of our three full 
geometry models (from Figure 2). The value obtained by Gulick was converted from 
volume to mass for the purpose of this figure, with an assumed water density of 1000 
kg.nf 3 . 

For magma chambers with aspect ratios greater than the baseline model, it is no 
longer reasonable to assume that heat is lost through the chamber walls alone, and the 
total discharges for our geometry and that of Gulick (1998) are expected to diverge. 
Instantaneous discharge diverges even at the baseline aspect ratio, with a maximum 
value about 1 .8 times higher for our geometry, and a steeper subsequent decay. 

Magma chamber cooling times are of interest in these models and they are 
defined for our finite difference model as the average time taken for chamber cells to 
cool below a specified threshold temperature. For all three models, the cooling time 
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for a threshold of 450 ° C (i.e. half the emplacement temperature) is approximately 30 
kyr. While the magma chamber is impermeable, the groundwater flow past its walls 
does not significantly enhance cooling (see Norton and Knight, 1 977). Once the 
chamber becomes permeable, however, cooling progresses at different rates in each 
model. The time taken for the chamber to cool to 250 °C is 67, 61, and 40 kyrs for 
permeabilities of 10' 17 , 10' 16 , and 10' 15 m 2 respectively. 



time [kyr] 


Figure 2.4 Surface discharge from magma chambers of volume 50, 100, and 2000 
km 3 , all at a depth of 8.5 km, an aspect ratio of 2, and with host rock permeability of 
10' 16 m 2 . The total discharge of the 50 km 3 model appears above its curve. The 
numbers above the other curves denote the fraction of this discharge that their 
corresponding models produced. 
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2.3.2 Volume 

Head and Wilson (1994) suggest that martian magma reservoir volumes could be 
as great as 2000 km J , and that chamber depths are most likely to range from 8 to 12 
km. we ran models of D/H = 2 magma chambers with volumes of 100 and 2000 km J 
at a depth of 8.5 km, in host rock with a permeability of 10‘ 16 m 2 (Figure 2.4). The 
100 knT chamber produced about 2.5 times as much discharge as the 50 km J 
chamber, while the jump from 50 to 2000 km J resulted in a factor of 1 80 increase. 
This approximately linear relationship is reflected in the magnitude of the 
instantaneous discharge, whose peak has proportional increases. 

The total (cumulative) discharge from the system scales with three independent 
factors: chamber volume, surface area-volume ratio A/V (1.59, 1.26, 0.46, with 
increasing volume), and the average background host rock temperature into which the 
chamber is emplaced. Chamber volume is a measure of the total energy available to 
drive hydrothermal circulation; A/V is a measure of the relative duration of cooling 
(and varies with volume despite the conservation of aspect ratio); and average 
background host rock temperature is a measure of small differences in cooling time 
due to the low geothermal gradient. The deeper the chamber floor, and therefore the 
higher the average ambient temperature, the smaller the initial temperature gradient 
across the chamber surface. This slows the flux of heat into the country rock, 
producing greater cooling times and discharge longevity. Volume, A/V, and 
temperature together predict cumulative discharge s (relative to the baseline model) 
of 2.5 and 150 for volumes of 100 and 2000 km 3 respectively, close to the modeled 
values. 

Although not readily apparent on the semi-log plot of Figure 2.4, the discharge 
produced by the 2000 km J chamber drops off at about 0.8 Myr, in agreement with 



values obtained by Cathles et al. ( 1 997) for a 2500 km 3 chamber with host rock 
permeabilities ranging from 4 x 10' i7 to 10' 16 m 2 . 
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Figure 2.5 Surface discharge from 50 km 3 magma chambers at depths of 2, 8.5, and 
15 km, all with aspect ratio 2, and with host rock permeability of 10' 16 m 2 . The total 
discharge of the 2 km deep chamber appears above its curve. The numbers above the 
other curves denote the fraction of this discharge that their corresponding models 
produced. 


2.3.3 Depth 

Models with 50 km 3 magma chambers at depths of 8.5 and 1 5 km and with a host 
rock permeability of 10' 16 m 2 were run (Figure 2.5). Total discharge (relative to the 
baseline value) is 0.69 and 0.42 for chamber depths of 8.5 an 15 km respectively. 
Similar depth dependence exists for the same three chamber depths at host rock 
permeabilities of 10' 17 and I0' 13 m 2 . 
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An approximate prediction of the change in total discharge can be made by 
accounting for changes in chamber depth and the average background host rock 
temperature into which the chamber is emplaced. The depth of the chamber floor 
(which therefore also accounts for chamber height) and the average background host 
rock temperature together predict total discharge (relative to the baseline value) of 
0.52 and 0.40 for chamber roof depths of 8.5 and 1 5 km repsectively, similar to 
modeled values. 

All three chambers cool to half their emplacement temperature in approximately 
30 kyr and further cooling is required to demonstrate the effect of the small 
geothermal gradient: the 2 km deep chamber cools to 250 °C in 61 kyr, and the 15 
km chamber takes 71 kyr. A summary of the total discharge of all nine depth and 
permeability combinations is depicted in Figure 2.6. 



log total 
discharge [kg] 


Figure 2.6 Total discharge from 50 km J magma chambers of aspect ratio 2, at depths 
of 2, 8.5, and 15 km, and host rock permeabilities of 10' 17 , 10" 16 , and 10' 15 m 2 . 
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2.3.4 Aspect Ratio 

Magma chambers with D/H = 0.2, 2, and 20 were modeled (Figure 2.7). These 
correspond to chamber radii of 1.2, 2.5, and 5.4 km, and chamber heights of 12, 2.5, 
and 0.54 km respectively. Chamber depth in all cases is 2 km, and chamber volume 
is 50 km 3 . The radial extent of the host rock was fixed at 50 km for all three models 
in order to observe and compare the scale of flow produced by each magma chamber. 
In all three models the flow pattern consists of a single, clockwise rotating convection 
cell alongside the chamber. For D/H = 20, this pattern does not give way to a series 
of cells above the chamber roof, as one might expect for flow between two infinite 
horizontal surfaces at different temperatures. 



time [kyr] 


Figure 2.7 Surface discharge from 50 km 3 magma chambers with D/H = 0.2, 2, and 
20, all at a depth of 2 km, and with host rock permeability of 10' 16 m 2 . The total 
discharge at D/H = 2 appears above its curve. The numbers above the other curves 
denote the fraction of this discharge that their corresponding models produced. 
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The discharge of these models is controlled primarily by the cooling time of the 
magma chamber and therefore its surface area-to- volume ratio A/V (1.89, 1.59, and 
4.06 km' 1 with increasing aspect ratio). Ratios of A/V predict cumulative discharges 
(relative to the baseline value) of 0.84 and 0.39 for aspect ratios of 0.2 and 20, again 
close to modeled values. The effects of aspect ratio are complicated by the presence 
of a substantial geothermal gradient, as discussed below. 

2.3.5 Geothermal Gradient 

Models with magma chamber aspect ratios of 0.2, 2, and 20 were run with an 
initial host rock geothermal gradient of 20 °C/km, perhaps representative of early 
Mars (Schubert et al., 1992). The flat, sill-like chamber (D/H = 20) produces the 
greatest peak discharge because of its large horizontal exposure. The position of this 
chamber in relatively cool, high elevation host rock results in a greater temperature 
gradient across its surface and therefore more rapid cooling. Conversely, the tall, 
pipe-like chamber (D/H = 0.2), being immersed in warmer temperatures at depth, 
cools more gradually. Its geometry produces the largest convection cell and offers 
minimum obstruction to flow, resulting in the greatest total discharge despite its low 
peak value. Overall, tenfold variations in aspect ratio produce changes in discharge 
of less than a factor of 3. 

The times taken for magma chambers to cool to half their emplacement 
temperature are also affected by geothermal gradient and are, in order of increasing 
aspect ratio, 25, 45, and 4.0 kyr respectively. (In the absence of a significant 
geothermal gradient, the same chambers cool in 16, 28, and 3.5 kyr respectively.) An 
8.5 km deep chamber cools to half its initial temperature in 52 kyr (as opposed to 30 
kyr with a negligible geotherm), and produces about 7 times as much discharge with a 
peak 3.5 times higher. 
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It should be noted that for permeabilities greater than 10' 16 m 2 , the Rayleigh 
number of a plane porous medium (e.g., Turcotte and Schubert, 1982) at 20 °C/km 
indicates that weak “free” convection may occur in the absence of a magma chamber. 
This phenomenon is indeed observed in high geothermal gradient models with no 
magma chamber. It has also been invoked by Raffensperger and Garven (1995) to 
explain the location of uranium ore deposits in sedimentary basins in Canada and 
Australia, and Travis et al. (2001) showed that it may be capable of melting 
significant volumes of subsurface ice in the martian crust. However, all of these 
models, including our own, have homogeneous permeability in the convecting zone: 
realistic vertical heterogeneity in planetary crusts will probably inhibit the 
development of large-scale free convection. While background geothermal gradients 
may have some affect on intrusion cooling times, we view their contribution to 
hydrothermal circulation as marginal. 


2.3.6 Ice 

Gulick (1998) suggested that a 2 km thick subsurface permafrost layer above a 50 
krrF magma chamber would melt in a few tens of thousands of years. This is of the 
same order as the time taken for discharge to peak in our own models, indicating that 
ice could significantly reduce total discharge, we thus ran the baseline model with a 1 
km thick layer of subsurface ice (Figure 2.8), modeled after Bonacina et al. (1973) 
who approximated melting by augmenting the specific heat over a small temperature 
interval AT. The augmmented value is given by 

L 
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where subscripts S and L refer to the solid and liquid phases respectively, and L is the 
latent heat of fusion, we used c ; . =1000 J.kg'.K' 1 , L = 3.34x 10 5 J.kg'\ and T=5°C, 


giving q. = 6.78x1 0 4 J.kg'.K' 1 . 
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Figure 2.8 Baseline model with an initially 1 km thick subsurface ice layer shown at 
15 kyr (A) and 59 kyr (B). The scale in each image is about 5 by 5 km. Water and 
steam velocities are shown by black and gray streamers respectively. Only every 
other row and column are shown. Shaded regions (apart from the ice) represent 
different fluid phases as described in the key. At 59 kyr, only a small remnant of 
steam remains and the magma chamber has become as permeable as the host rock. 


The hydrothermal system took 52 kyr to melt a hole in the ice. The total discharge 
was 2.95 x 10 11 kg, about a quarter that of the baseline model. This reduction (which 
is expected to be worse for deeper chambers and lower host rock permeabilities) may 
be useful for constraining an upper bound on permafrost thickness during valley 
formation on Mars, if indeed any permafrost is present during the Noachian. For 
valleys to form through sapping processes alone, the permafrost must be melted 
through, or aquifers carrying groundwater beneath the permafrost must intersect the 
surface. In either case, the permafrost can be no thicker than a few hundred meters. 
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2.3. 7 Water Table 

Until now we have assumed a completely saturated host rock. If an unsaturated 
zone is present, part of the energy in the system must go toward elevating water to the 
surface, with the result that discharge is diminished, we compare with hydrostatic 
conditions the vertically integrated product of density, gravitational acceleration and 
depth in order to estimate the height a water table may attain through thermal- 
convective expansion. Integrations are performed over all vertical columns of finite 
difference blocks at all times, and for a range of water table depths. The discharge is 
calculated from those parts of the water table that intersected the surface, and 
compared to zero depth water table models. The results (Figure 2.9) show that 
models with a water table 50 m deep produce only 10% of their zero depth water 
table equivalents, while models with a 100 m deep table produce only 2%. 
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Water Table Depth [m] 


Figure 2.9 Discharge as a function of water table depth, normalized to the zero depth 
water table value. 


2.4 Discussion 

2.4. 1 Implications for Hydrothermal Alteration 
Water/rock ratio (W/R) can be an important influence on the mineralogy of 
alteration products in hydrothermal systems. Martian hydrothermal alteration is 
thought to occur at low W/R (< 10 by mass) (Griffith and Shock, 1997, and Newsom 
et al., 1999), in which case its effect is small, and the initial mineral composition is 
the primary influence on the alteration assemblage (Griffith and Shock, 1997). 

HYDROTHERM mass flux results allow the W/R of a hydrothermal model to be 
calculated. The total flux of water through regions of the model above a specific 
threshold temperature are integrated over time. In this way, the W/R, or the total 
water mass that makes contact with the reacting rock, can be calculated for individual 
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finite difference cells. The W/R value averaged over all cells is a suitable 
representative of the system W/R, and is shown for the three reaction temperatures 
modeled by Griffith and Shock (1997), and for various model dimensions, in 
Table 2.1. 


Threshold 
temperature ( ° C) 

k = 10' 15 rrf 
volume = 50 km 3 

k = 10‘ 16 m 2 
volume = 50 km 3 

k = 10' 16 rrf 
volume = 2000 km 3 

150 

0.12 

0.025 

0.067 

200 

0.18 

0.037 

0.10 

250 

0.28 

0.054 

0.15 


Table 2.1 Water/rock ratios (by mass) calculated for 8.5 km deep chambers of the 
indicated volumes and host rock permeabilities k. 


The large velocities in the k= 10' 15 m 2 model give the highest W/R for all 
threshold temperatures, approximately 5 times that of the A: = 10' 16 m 2 model with a 
chamber of the same volume. Negligible geothermal gradient and low host rock 
permeability in our models are both factors that contribute to small W/R. Typical 
martian geotherms (~20 °C/km) will lead to larger W/R, so these estimates can be 
viewed as lower bounds. The presence of a significant geotherm may be expected to 
increase discharge by a factor of about 3 or 4 (as observed in the model with an 8.5 
km deep chamber with geothermal gradient), and since W/R is directly proportional 
to the discharge flowing through the alteration area, it may experience a similar 
increase. This would still leave most values listed in Table 2.1 below unity. 

A primary goal of the work by Griffith and Shock (1997) was to model the 
amount of water bound to the rock during alteration. In a model based on the 
Shergotty SNC meteorite, about 8% of the final mineral assemblage (by weight) was 
water. It follows that in reactions with a W/R of 0.08, all of the water coming into 
contact with the rock would be sequestered in the final mineral assemblage - 
additional water for surface discharge would require W/R’s in excess of 0.08. 
Assuming the Shergotty composition is sufficiently generic to be applied to our own 
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models, a W/R limit of 0.08 implies a minimum permeability of 10' 16 nr fora 150 °C 
reaction in host rock surrounding a 2000 km 3 chamber (Table 2.1). A caveat to this 
conclusion is that some water driven to the surface by hydrothermal convection will 
not pass through zones of chemical alteration, and will thus avoid being sequestered 
in the host rock. 

Newsom et al. ( 1 999) suggest that the relative abundances of mobile elements 
such as sulfur, chlorine, sodium, and potassium may be explained by the presence of 
a mixture of neutral-chloride and acid-sulfate fluids during soil formation. 

Production of the latter fluid requires vapor transport (Rye et al ., 1992), and 
Ingebritsen and Sorey (1988) discuss situations in which vapor-dominated zones may 
occur. Their models require combinations of low-permeability barriers and in some 
instances topographic gradients to sustain vapor-dominated zones, which develop in 
the shallow subsurface only. These specialized structures have not been included in 
our generalized models; nonetheless, our highest permeability model (10~ 15 m 2 ) with 
the shallowest chamber (2 km) does produce a short-lived (few thousand years) two- 
phase zone between magma chamber and surface. Steam develops here because of a 
combination of low pressures (which drop to a minimum of about two thirds 
hydrostatic pressure) and high temperatures. At lower permeabilities, fluid rising 
from the chamber does not transport enough heat upwards to elevate temperatures to 
the required levels for steam production. 

If substantial quantities of vapor are required for soil alteration, then the 
abundance of steam in our high-permeability model may point to 10' 15 m 2 as a lower 
bound to permeability, a somewhat tighter constraint than that imposed by W/R ratio 
alone. Alternatively, if volcanic aerosols (Newsom et al., 1999) or basalt 
palagonitization (McSween and Keil, 2000)) produce sufficient soil alteration, then 
no such limit is necessary. 
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2.4.2 Relationship of Hydrothermal Circulation to Valley Networks and Magnetic 

Anomalies 

The apparent spatial correlation between valley networks and magnetic anomalies 
suggests a link between processes involving the acquisition of thermoremnance in a 
cooling intrusion and processes involving the production of surface water for the 
creation of fluvial channels, i.e. hydrothermal circulation. An overlay of maps 
(Figure 2.10) showing the valley networks (Kieffer, 1981) and the vertical component 
of the magnetic anomalies as measured by the Mars Global Surveyor magnetometer 
(Acuna et al., 1999) inferred at a constant altitude of about 200 km (Purucker et al., 
2000) visually suggests that the two are correlated (Jakosky and Phillips, 2001). 



150 



nT 


Figure 2.10 Overlay of the absolute values of the 200 km vertical magnetic 
anomalies (Purucker et al., 2000) and the valley networks (Kieffer, 1981) (cylindrical 
projection). Valleys were not mapped for latitudes below 60°S. The anomalies range 
from 0 to 670 nT. 
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A statistical analysis of the relative positions of valley networks and magnetic 
anomalies quantifies the correlation. A binomial test was performed, with success 
defined as the occurrence of a valley network and a magnetic anomaly (with vertical 
component above a specified threshold) in the same data bin. we consider only the 
southern highlands in these calculations. The valley network map (Kieffer, 1981 ) 
contains the coordinates of the 0.25 by 0.25 degree bins that contain valleys; the 
number and length of the valleys are not quantified. The magnetic anomaly map 
(Purucker et al., 2000) is binned at 1 by 1 degree and was regridded at 0.25 degrees to 
match the valley grid. The number of degrees of freedom must, however, be adjusted 
to match the true spatial resolution, which is limited by the magnetic data. The total 

number of 0.25 degree bins is therefore multiplied by the quantity 


e ~ 


'0.25V 

V ^ m J 


where d m is the length scale in degrees for magnetic resolution, we determined that 
the variogram for the 200 km magnetic field reaches half of its sill (asymptotic) value 
in 200 km and 90% of the sill value in 400 km, in agreement with the rule of thumb 
for potential fields that the spatial resolution is approximately equal to measurement 
altitude. Therefore d m ~ 200 km ( « 3 degrees) may be most appropriate, but we 
consider a range from 2 to 6 degrees. 

Relevant probabilities may be obtained using the total number N of 0.25 by 0.25 
degree bins in the region of interest, the number n of these bins containing valleys, 
the number m of bins with magnetic anomalies above a specified threshold, and the 
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number C of correlations. The observed correlation probability is then p-C In, 
while the expected correlation probability for a valley placed randomly in the region 
of interest is po = m / N. The probability q of obtaining C or more correlations if the 
n valleys are placed randomly in the region of interest, can be obtained using the 
cumulative binomial distribution with sC, en, and poas the minimum number of 
successes, the number of independent trials, and the probability of a single success, 
respectively. A low value of q implies that the relative distribution of valleys and 
magnetic anomalies on Mars is not what one would expect from chance. 

The results for a range of magnetic anomaly thresholds (Table 2.2) indicate that 
minimum chance probabilities occur near a 10 nT threshold for the vertical magnetic 
field at 200 km altitude. The area containing magnetic anomalies falls off sharply at 
higher thresholds, allowing a greater probability of chance correlation. At lower 
thresholds, the entire map is considered anomalous and so again it is easier to produce 
randomly the observed correlation. Over length scales at which the magnetic 
anomalies can be described as strongly coherent ( d,„ < 4 degrees), the probability of a 
chance correlation is relatively small (q < 0. 16). 


magnetic 

threshold 

[nT] 

Po 

P 

(P~Po)fp d m = 2° 

<4 = 3 0 

<7 

<4 = 4 0 

<4 = 6° 

1 

0.896 

0.962 

0.069 

0.0785 

0.182 

0.494 

0.464 

2 

0.827 

0.928 

0.109 

0.0204 

0.102 

0.210 

0.654 

5 

0.669 

0.827 

0.191 

0.00228 

0.0480 

0.172 

0.268 

10 

0.469 

0.658 

0.287 

0.00223 

0.0345 

0.160 

0.434 

20 

0.295 

0.465 

0.366 

0.00459 

0.0590 

0.166 

0.343 

50 

0.126 

0.201 

0.373 

0.0686 

0.303 

0.330 

0.612 

100 

0.0521 

0.0818 

0.363 

0.260 

0.450 

0.575 

1 


Table 2.2 Statistical data for the correlation between valley networks and magnetic 
anomalies. 


A genetic correlation requires that the valley networks and magnetic anomalies 
are the same age. A significant majority of valleys are Noachian (70-92%; Scott and 
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Dohm, 1 992; Carr, 1 996); many of those that are younger are not in the southern 
highlands and so are excluded from our survey. The magnetic anomalies, because of 
their inferred deep-crustal origin (see Proposition 1 below) without surface 
manifestation, must also be very old. However, even the oldest valley networks 
individually preserve only some part of the Noachian that was not subsequently 
locally erased, whereas the magnetic anomalies probably reflect a greater span of 
crustal history. In other words, some valley networks that were associated with 
magnetic anomalies may have been resurfaced, whereas other valley networks may 
have formed subsequent to emplacement of the magnetic anomalies by hydrothermal 
or other processes. The normalized excess of correlated valleys and magnetic 
anomalies (p-po)/po (Table 2.2) may be taken as representative of the fraction of 
valleys for which a genetic correlation may be inferred, between about one-quarter 
and one-third. 

The role of hydrothermal circulation in the relationship between valley networks 
and magnetic anomalies may now be described in terms of a central hypothesis 
composed of two main propositions, defined and discussed in the following sections. 

Proposition 1. The first proposition is that the magnetic anomalies formed as 
intruded crust and that the acquisition of thermal remnant magnetization (TRM) 
occurred at relatively great depth. The strong observed magnetizations in the martian 
crust of 20 to 40 A.m 2 (Connemey et al., 1999; Grimm, 2000) imply magnetization 
depths of up to a few tens of kilometers. The presence of magnetic anomalies in 
Arabia Terra, which has been strongly resurfaced (McGill, 2000; Hynek and Phillips, 
2001) also points to a deep origin. On Earth, Layer 3 gabbros of the oceanic crust are 
magnetized to the extent that they contribute between 25% and 75% of the observed 
marine anomalies (Pariso and Johnson, 1993), while the underlying mantle is 
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unmagnetized (Wasilewski et al., 1979). Therefore deep crustal magnetization of 
Mars is reasonable. 

The mineral composition of the magnetized material producing the magnetic 
anomalies is likely to contain magnetite or hematite as the primary magnetic carrier. 
Although magnetite is generally favored, there is much support for hematite (e.g., 
Connemey et al., 1999). Kletetschka et al. (2000) show that for an applied magnetic 
field of 0.1 mT (about twice the strength of the earth’s present geomagnetic field), 
multidomain hematite reaches maximum TRM saturation, whereas magnetite only 
reaches only a few percent thereof. 

Proposition 2. The second proposition is that hydrothermal discharge attending 
crustal formation processes in the southern hemisphere of Mars was sufficient to 
provide the water necessary to carve the valley networks. Hydrothermal systems on 
earth exhibit such discretized concentration of outflow, as evidenced by the presence 
of geysers and springs rather than diffuse outflow everywhere above the magmatic 
intrusion. The martian valley networks are characterized by low drainage densities, 
implying again that crustal heterogeneities may localize discharge. 

Testing this proposition requires knowledge of the amount of water necessary to 
erode the valley networks, the total amount of water available to hydrothermal 
systems, and the actual hydrothermal discharge produced. An estimate of the 
required water volume for valley erosion may be made from values of areal coverage, 
drainage density, valley cross-section, and sediment-to- water ratio. Using a map of 
drainage densities (Carr and Chuang, 1997) we estimate the total area covered by 
valley networks to be about 1 .4 x 10 7 km 2 . Since we estimated earlier that only one 
quarter of the valley networks may preserve direct interaction with hydrothermal 

i ■ 'y 

systems, we use a reduced effective area of 3.6 x 10 km'. Carr and Chuang ( 1 997) 
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calculated a globally averaged drainage density of 0.0032 km' 1 . The product of 
effective area and drainage density, multiplied by typical valley width and height (5 
km and 150 m respectively), result in 8.6 x 10 J km J of removed material. Sediment- 
to-water ratios may range from 1:4 to 1 : 1000 (Gulick, 1998, and references therein), 
implying that the volume of water required to erode the valleys was between 3.5 x 
10 4 and 8.6 x 10 6 knT, equivalent to a global water layer between 0.2 and 60 m deep. 
These values are well below the estimate of hundreds of meters for the global crustal 
inventory, indicating that discharge from the valleys had a negligible impact on the 
global water budget. They further imply that discharged water need not have been 
recharged to the crust. 

The total hydrothermal discharge produced is computed as follows. First, we 
assume that each magma chamber contributing to crustal formation intrudes into 
steady ambient temperature and pressure conditions; this is most likely if intrusions 
that formed the southern highlands moved between different loci rather than 
spreading from a single location (Grimm, 2000). In this way, the discharge 
contribution from a single intrusion is just that of its equivalent HYDROTHERM 
model, we assume further that the crust covering the area occupied by valley 
networks (as calculated above) is eventually built up to a depth of 20 km by magma 
chambers packed side by side and one on top of the other, we sum the total 
discharges from each intrusion to calculate the total mass of surface water produced, 
we consider only those intrusions that contribute to magnetic anomalies, discarding 
other intrusive events. Determining the relative contribution of individual intrusions 
is not possible, but the probability po may be used as an indicator of the appropriate 
fraction to be considered. Its value for a magnetic threshold of 10 nT, i.e. 0.469 
(Table 2.2), is applied to all of our results. 
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Figure 2.11 Summary of global discharge production (by mass, in kg) for various 
model parameters. The two thick dashed lines indicate the bounds on the mass of 
water required to erode the valley networks. Represented are discharge values for 
different magma chamber aspect ratios, host rock permeabilities, magma chamber 
volumes, and for the baseline model with different assumed water table depths. Note 
that the baseline model (D/H = 2 ,k = 10' 16 m 2 , chamber volume = 50 km J , water table 
depth = 0 m) appears in each vertical group of bars. 


A summary of discharge production over the entire region of interest is shown for 
different models in Figure 2.1 1. The two thick dashed lines denote the bounds on 
total required discharge. Represented are discharge values for different magma 
chamber aspect ratios (all with a volume of 50 km J , host rock permeability of 10' 16 
m 2 ), host rock permeabilities (all chambers have D/H = 2 and a volume of 50 km 3 ), 
magma chamber volumes (all chambers have D/H = 2; host rock permeability is 10' 16 
m 2 ), and values for the baseline model with different assumed water table depths. 
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All models meet or exceed the discharge production requirements, although it 
should be noted that if the k = 10' 17 m 2 model is assumed to have a water table of 100 
m, its production will fall below the minimum required value. Other factors such as 
evaporation may further reduce the discharge available to carve valleys. The present 
evaporation mass flux is likely in the region of 4 x 10' 8 kg.m' 2 .s'' (Ingersoll, 1974), 
which would reduce the effective discharge by as much as an order of magnitude or 
more. Ice may also reduce discharge (Section 3.6). A 1 km thick layer of subsurface 
ice in our baseline model causes a 75% drop in total discharge. In models with other 
parameter values (e.g. greater magma chamber depth and smaller host rock 
permeability), ice may inhibit discharge more severely, if not completely. 

2.5 Conclusions 

In numerical models of martian hydrothermal systems, we explored the control on 
surface discharge of magma chamber depth, volume, aspect ratio, and host rock 
permeability and porosity. Discharge has an approximately linear relationship to 
magma chamber volume and host rock permeability (in the range 10" 17 to 10' 15 nr). 
The influences of depth and aspect ratio are weaker, and that of porosity is negligible. 

Some geochemical aspects of martian hydrothermal systems were considered by 
calculating water/rock ratios in our numerical models at various reaction 
temperatures. Ratios tend to be low, but sufficiently large at mean permeabilities > 
10‘ 16 nr for groundwater flow to be sustained, consistent with the expected storage of 
water in alteration assemblages. The presence of a short-lived vapor-dominated zone 
in our model with high host rock permeability (10' 15 m 2 ) and shallow chamber depth 
(2 km) suggests that hydrothermal alteration processes may be responsible for the 
observed relative abundances of certain salts in the martian soil, although other forms 
of alteration should not be excluded. 
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Crustal formation processes which formed the magnetic anomalies observed on 
Mars today may have been attended by hydrothermal circulation that also provided 
surface water for valley network erosion. This idea is in agreement with the observed 
spatial correlation between magnetic anomalies-and valleys and was tested further 
within the framework of a central hypothesis made up of two propositions. The first 
is that the magnetic anomalies formed as intruded crust and that the acquisition of 
thermoremnance occurred at relatively great depth. The second is that hydrothermal 
discharge attending global crustal formation processes is sufficient to provide the 
water necessary to carve the planet's valley networks. 

we tested this hypothesis using the numerical models described above, assuming 
that the crust was formed by the heterogeneously spaced and timed intrusion of 
multiple magma chambers, each of which produced the discharge predicted by its 
individual numerical model, we determined that many model configurations in the 
explored portion of the parameter space were capable of producing sufficient water to 
erode those valley networks statistically related to hydrothermal circulation. In 
particular, modest crustal permeabilities of 10' 16 to 10' 15 m 2 can produce the discharge 
required to carve valleys and satisfy geochemical constraints, even in the presence of 
mitigating factors such as evaporation and finite water-table depth. 
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CHAPTER 3 

LONG RUNOUT LANDSLIDES 

3.1 Introduction 

Landslide features on Mars bear an unambiguous resemblance to terrestrial 
landslides. They occur mainly in the equatorial Valles Marineris canyon system, in 
both open and closed depressions (Lucchitta 1979). The aureole deposits on the lower 
flanks of the largest Martian shield volcano, Olympus Mons, may also be landslides 
(Lopes 1980, 1982), although slower processes such as gravity spreading may have 
been responsible for their formation (Francis and Wadge 1983, Tanaka 1985). In 
many Valles Marineris landslides a curved scar cut into the canyon rim indicates the 
site of the original failure. The scar is typically a few tens of kilometers wide and 
recedes a few kilometers into the canyon wall. Landslides occur in walls as high as 8 
km and cover horizontal distances (runouts) as great as 80 km if unimpeded by 
obstacles. The volume of individual landslide deposits ranges from tens to thousands 
of cubic kilometers. 

Of ail the processes that have shaped the Valles Marineris, landsliding is among 
the most recent. The landslide scars are smooth and the deposits are sparsely cratered 
and show little evidence of deflation (Lucchitta 1992) suggesting a late Hesperian to 
Amazonian genesis. Slide scars cut into pre-existing spur-and-gully formations, and 
the slide material overlies interior deposits. There is, however, some indication that 
tectonic movement postdating landslide deposition has occurred along faults closely 
related to the formation of the Valles Marineris (Lucchitta 1979, Lucchitta et al. 

1994). Earlier marsquakes triggered by tectonic activity may have been responsible 
for triggering the landslides (Caruso, 2003). 
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Although reentrants are generally smooth, their upper elevations are typified by 
dark vertical banding which, on close inspection of MOC images, is seen to consist of 
small, rugged spur formations, sometimes horizontally striated, that are consistent 
with the volcanic nature of the plateaus surrounding the Valles Marineris (Scott and 
Tanaka 1986). Emanating from the lower boundaries of these miniature spurs are 
lighter colored, linear, overlapping talus deposits blanketing the lower elevations of 
the scar. The influence of material strength heterogeneity on landslide dynamics is 
discussed below. 

The landslide deposits themselves have small-scale features, some of which are 
difficult to explain. These include longitudinal grooves, rather than transverse ridges 
like those observed on terrestrial slides. Lucchitta (1979) suggested that longitudinal 
grooves result from differential forward velocities that can be explained by the 
presence of a fluid. Many large deposits, both martian and terrestrial, have runouts 
that imply substantially lower internal friction angles than those typical of rocky 
materials. Processes thought to produce long runouts typically require fluidization by 
water (Habib 1975, Goguel 1978) or air (Kent 1966, Shreve 1968). Arguments have 
been given in support of both the presence (Lucchitta 1987) and absence (McEwen 
1989) of water in the Valles Marineris slides. Chaotic terrains at canyon termini, 
fluvial features on canyon rims (Carr, 1996; Coleman et al., 2003), and the possiblity 
of ice-related processes in interior layered deposit formation (Weitz and Parker, 2000) 
all point toward the past presence of water in the Valles Marineris. 

The thin Martian atmosphere is unlikely to be dense enough to provide 
lubrication, but outgassing of carbon dioxide from the soil on failure could play an 
important role. Some fluidization processes do not require an actual fluid: in 
dispersive grain flow (Bagnold 1954), the weight of the material is supported by 
impacts between its constituent rock particles, whereas for acoustic fluidization 
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(Melosh 1979) an acoustic field transmitted through the material is sufficient to 
support for brief periods the static overburden pressure without requiring separation 
of particles. 

Features of landslide morphologies may be used to infer their dynamic behavior. 
Rheology is an important influence on the final configuration of a landslide: given a 
particular runout path and initial mass configuration, it is largely the rheology that 
determines key features such as runout length and final deposit thickness and slope. 
Here we report the results of numerical models of martian landslides designed to 
constrain rheological parameters. Comparisons are made with terrestrial and lunar 
models. 


3.2 Model 

3.2.1 Numerical Simulations 

we use a dynamic analysis model (DAN, Hungr 1995) which, given the runout 
path and initial mass profile, computes the time-varying shape and velocity 
distribution of the slide. DAN has been successful at modeling terrestrial landslides 
with a wide range of rheologies, including sand and coal failures (Hungr 1995) and, 
where physical parameters are not available for comparison, it successfully 
reproduces other numerical analyses (Hungr 1995), including that of the Madison 
Canyon, Montana rock avalanche (Trunk et al 1986). 

A finite difference Lagrangian approach is used to track individual mass blocks as 
they move together down the runout path. Each slide is modeled as a single row of 
blocks (20 in our models) whose initial heights depend on the prescribed initial shape 
of the failed mass. Quantities that vary normal to the runout path are represented by 
their average values. For each block, the component of the gravitational force along 
(tangential to) the runout path works against a pressure contribution and a basal 
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resistance force. The pressure contribution arises from the differing heights of 
adjacent blocks. Material tends to shift laterally to reduce high, unstable surface 
slopes. The resulting communication of pressures throughout the landslide is 
controlled by a lateral pressure coefficient whose lower and upper bounds reflect the 
tensional and compressional strengths of the material respectively. The pressure 
coefficient has a small influence on the final profile and is not considered to be of 
primary importance. 

The basal resistance force for each block is determined by one of several available 
rheologies, which include frictional, Bingham, and power law. In the frictional 
rheology, the basal resistance offered by a block is proportional to the normal stress it 
exerts on the runout path. The constant of proportionality is the friction coefficient, 
often expressed as the tangent of the angle of internal friction <j>, set to 20° in all 
models unless where otherwise specified. By default, DAN uses the same static and 
dynamic friction angles. In reality, the dynamic friction angle is less than the static 
value (the latter is typically around 30° for fractured crustal rock, Brace and Kohlstedt 
1980) and may change with the velocity of the material (Jaeger and Cook 1979). 

Other workers (e.g., Wang 1 997) have used models with dynamic friction coefficients 
as low as 0.6 times the static value (i.e. a static angle of 30° corresponds to a dynamic 
value of 19°). There is also a dependence on clay content (Skempton 1964) that could 
reduce the dynamic angle by a few degrees. Small quantities of clay minerals have 
been observed indirectly in spectra (e.g., Soderblom 1992) and directly in the SNC 
meteorites (e.g., Gooding et al. 1991 and Treiman et al. 1993). Our choice of friction 
angle thus represents a physically plausible lower limit. The frictional law used in 
DAN is such that a change in friction angle will have the same affect on dynamics as 
a corresponding change in pore fluid pressure. A lower limit for friction angle thus 
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precludes the interpretation of pore pressure as artificial compensation for an 
exaggerated friction angle. 

The basal resistance also depends on the pore-pressure coefficient r u , defined as 
the ratio of pore pressure exerted by an interstitial fluid to the total normal stress at 
the base of the block. Pore fluid carries some of the weight of the overlying material, 
reducing basal friction, we stress that a given pore pressure coefficient describes 
conditions at the failure surface only and heterogeneities in permeability may allow 
other parts of the failed wall section to remain dry. Consequently, pore pressure 
coefficients that correspond to hydrostatic or even superhydrostatic pressures do not 
necessarily imply total saturation at all depths. The basal resistance is also affected by 
centrifugal forces, since these increase or decrease the effective normal stress 
according to the curvature of the path. 

The Bingham rheology, used successfully to model rock avalanches (Hungr 1995, 
and references therein) consists of a Newtonian fluid with a finite yield strength to. 
Assuming that shear stress increases linearly with depth, the constitutive equation 
(which relates shear stress to velocity gradient) may be integrated to yield an 
expression for the basal resistance in terms of fluid velocity v, yield strength, and 
dynamic viscosity p. (The basal resistance is simply the shear stress at the base of the 
slide.) The yield strength gives rise to a solid cap which rides on top of the flow. As 
the deposit extends itself along the runout path, its tensile strength slows it down, 
causing the cap to thicken until it consumes the entire flow, bringing the deposit to a 
halt. This effect is exacerbated by the lowering of the cap as the deposit thins. 
Landslides with higher yield strengths have thicker caps and so form thicker deposits 
with shorter runouts. 

A power law may be used to model flows in which most of the lateral shear is 
concentrated in a thin basal layer (Melosh 1987) such as the Blackhawk landslide in 
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California (Shreve 1968). The Elm landslide in Switzerland was simulated 
successfully by Hsu (1975) in a scaled down model using a bentonite suspension, a 
fluid with power law behavior (Besq et al 2000). For a power law fluid, the shear 
stress is 


dvT 
dz ) 


3.1) 


where z is depth and n is a unitless index. In Equation 3.1, the symbol p does not 
have units of real viscosity and is called the “apparent” viscosity, we adopt n = 0.125 
(suggested by Melosh (1987) for acoustic fluidization) for all slides and for both 
fluidization processes. A value of 1 represents a Bingham fluid, while a value of 0 
corresponds to a rigid sliding block, we investigated the sensitivity of runout to n for 
a martian slide (Slide 3 below) and found that within the range of 0 to 1, halving or 
doubling the value of n increases or decreases runout by 2 to 2.5 km for a slide with a 
runout of 72 km at n = 0.125. The uncertainty in our choice of n is thus largely 
immaterial. 

we consider next a general fluidization process in which the slide is frictional on 
failure but becomes fluidized when sufficiently energetic, behaving more like a power 
law fluid. Lucchitta et al. (1992) explicitly suggests such a rheological transition. As 
noted above, pore pressure may be specified explicitly in the frictional rheology. With 
the power law, however, water content is implicit and its influence is assumed to 
manifest itself through the determination of yield strength and apparent viscosity. 

A special case of general fluidization is acoustic fluidization (Melosh 1979), 
which can occur without the presence of a liquid or gas and whose frictional stage 
may thus be associated with zero pore pressure. When acoustic waves of sufficient 
strength are transmitted elastically through the individual constituents of the landslide 
material, intra-particle rarefractions produce short periods of decreased overburden 
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pressure, temporarily reducing basal resistance. The acoustic field builds up to the 
required strength some time after the initial failure, so fluidization does not occur 
immediately. Although acoustic fluidization may be possible in the presence of an 
interstitial fluid, we address here the case of dry material only. This is, in fact, the 
fundamental distinguishing feature of this rheology and since it includes no explicit 
dependence on the acoustic processes described above, it may encompass other 
processes in which frictional and power law behaviours occur in the absence of fluid 
pressure, if others indeed exist. 

For clarity we note that the general and acoustic fluidization processes are not 
modal rheologies, but sequential combinations of two rheologies: frictional and 
power law. In general, the yield strength and apparent viscosity of the power law 
stage for a particular slide will not be the same for acoustic fluidization and general 
fluidization. The values of these parameters depend on the position at which the 
frictional rheology ends and the power law begins, and also on the velocity and stress 
distribution of the material at the transition. Transitions occurring further along the 
runout path have a frictional stage ending with different velocities, and the power law 
stage must thus have a different yield strength and/or viscosity in order to conserve 
runout. 

Inputs common to all rheologies include unit weight y, volume yield rate E (a 
measure of material deposition or entrainment, described below), and the lower and 
upper bounds of the lateral pressure coefficient. Unit weight was taken to be 10.00 
kN.m' 3 for martian slides, 26.44 kN.m' 3 for terrestrial slides, and 4.371 kN.m' 3 for 
lunar slides, corresponding in all cases to a density of 2700 kg.m' 3 . 
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3.2.2 Landslide Geometry 

Landslide deposit volume Vwas estimated for martian slides by summing 1/32 by 
1/32 degree Mars Global Surveyor (MGS) MOLA elevation data over the area 
covered by the material (Viking images were used to determine the location of 
deposit boundaries). For these calculations, the base of the deposit was assumed to be 
horizontal and its elevation was inferred from the surrounding canyon floor, or 
matched to the toe of the deposit. Our method of calculating volumes differs from 
previous approaches (e.g., McEwen 1989) which used the scar geometry to estimate 
the initial volume of the failed section of canyon wall. Since the configuration of the 
pre-slide wall is unknown there is considerable uncertainty in measuring volume this 
way; a direct measurement of the final deposit volume is likely to be more accurate. 

we assumed initially that volume was conserved during the slide event (i.e. no net 
deposition or entrainment of material). Approximate measurements of scar width W, 
scar slope ai, and typical nearby wall slope 0,2 allowed us to infer a suitable initial 
configuration for the failed material, with the required volume (Figure 3.1). Our use 
of volume as the primary measurement used to reconstruct the initial configuration is 
in line with experimental evidence (Hsu 1975) showing that volume, rather than 
geometry, has the stronger influence on runout. Nonetheless, the initial configuration 
was adjusted to match the shape of the unfailed wall adjacent to the landslide scar, 
where possible. Similarly, the runout path was shaped according to nearby wall and 
floor geometry. The coordinate system (Figure 3.1) was positioned such that the back 
edge of the top surface of the initial shape was at x = 0. The vertical coordinate 
represents elevations relative to the MOLA datum. 
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Figure 3.1 Conceptual model of initial landslide shape. Symbols correspond to those 
defined in the text. Runout occurs in the positive x direction. 

MOLA data were used to obtain vertical cross-sections of the landslide deposits 
along the direction of motion (all cross-sections mentioned henceforth refer to this 
plane). These were compared to DAN results to assess the suitability of the chosen 
rheological parameter values. Values were iteratively adjusted and the model rerun 
until a suitable final profile was achieved. 

In DAN, different sections of the runout path may be assigned different 
rheologies, we took advantage of this feature to simulate general and acoustic 
fluidization. The dynamics of the transition from frictional to power law is not well 
established and we chose the bottom of the initial steep section of the runout path as 
the default transition point. The transition position was varied about this point from 
model to model depending on the characteristics of the flow (such as velocity) and 
structural features observed in the final deposit profile (such as abrupt decreases in 
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slope). Rather than an abrupt change in rheology, this transition should be seen as the 
point at which the power law rheology becomes better than the frictional rheology at 
describing the motion. 

The conceptual model of landslides described thus far is somewhat limited by its 
essentially two-dimensional nature since changes in runout path width (dimension y 
in Figure 3.1) may influence the final cross-section of the deposit. In most cases, the 
landslide scar is curved in the cross-slope direction, i.e. each contour line is crescent- 
shaped. (In the downslope direction the scars, where exposed, are remarkably linear). 
This trend, which supposedly continues into the buried parts of the scars, suggests 
that material converges during the initial, steep part of its journey. On reaching the 
level, unconfined canyon floor, most of the energy in the slide has been channeled 
directly forward, resulting in a toe which travels further than the laterally spreading 
material on either side. The resulting deposit is shaped like an arrowhead pointing in 
the direction of motion. Despite lateral spreading, the average final deposit width is 
frequently smaller than the average initial width. Indeed, initial cross-sectional areas 
of the slides studied here are frequently smaller than final cross-sectional areas (Table 
3. 1 ), implying an average decrease in width in order to conserve volume. One can 
specify the required deposit width at each point along a DAN runout path. This 
feature is intended to simulate the effects of a laterally confining channel and is thus 
suitable for the curved landslide scar in our models. Unconfined spreading is best 
modeled by fixing the path width on the canyon floor to match the observed deposit. 

Another way to account for changes in cross-sectional area is to relax constraints 
on volume conservation. Landslides typically accumulate new material during the 
initial downhill phase. In DAN, one can specify the desired amount of entrainment or 
deposition over segments of the runout path. Individual block volume is updated at 
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each time step as follows (Hungrartd Evans 1997) 

V^V'j+EAs^- (3.2) 

where Fis the total slide volume, V) and Fj are the old and new volumes of the 
jth block, respectively, and As is the displacement of the y'th block during the time 
step. Parameter £ is a volume per unit runout path length and specifies how much 
material is deposited or entrained (depending on its sign). If the entire slide passes 
over a segment of the runout path with length a and entrainment rate E, the total 
change in volume is given approximately by aE. we assigned a finite entrainment rate 
to the initial, steep section of the runout path, and a value of zero elsewhere. Values 
of E for our models range from 10 4 to 10 7 m J .m'' depending on the initial volume of 
the slide and the required change therein. 


slide 

number 

H 

(km) 

H / L f 
(unitless) 

<Xi 

(°) 

D 

(km) 

initial cross- 
sectional 
area (km 2 ) 

final cross- 
sectional area 
(km 2 ) 

1 

1.1 

0.112 

19.0 

0.26* 

0.17 

0.15 

2 

8.0 

0.173 

24.2 

8.6 

74 

120 

3 

7.7 

0.0980 

24.2 

8.6 

67 

84 

4 

5.0 

0.0966 

14.1 

4.0 

27 

26 

5 

4.5 

0.122 

26.4 

0.47 

1.9 

2.2 

6 

5.6 

0.107 

16.4 

3.1 

16 

17 

7 

7.3 

0.101 

22.0 

7.4 

64 

105 

8 

10 

0.0145 

3.31 

86* 

720 

720 

9 

2.4 

0.273 

24.8 

0.022 

0.048 

0.048 


Table 3.1 Information inferred from vertical cross-sections parallel to the direction of 
motion. In some Valles Marineris slides, the initial mass was modeled after nearby 
unfailed canyon wall and the value for a 2 represents an average, rather than an exact, 
value. Depth values with an asterisk indicate the average value for a wedge-shaped 
initial shape. 


we ran models in which changes of cross-sectional area were produced 
exclusively by either width or volume variations, or a combination of both. In the last 
case, measurements of the landslide geometry were used where possible to constrain 
the width variations, and volume variations were adjusted to account for the 
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remaining change in cross-sectional area. The models presented here are exclusively 
of this last kind, although we note that models with only one kind of variation did not 
produce significantly different results. 


3,3 Results 

we present detailed models of a terrestrial landslide, six Valles Marineris 
landslides, part of the Olympus Mons aureole, and a lunar landslide. Four of the 
Valles Marineris slides have been described in detail by Lucchitta (1 979) and we 
refer to her landslide numbering scheme for easy reference. Images of Slides 1 to 9 
are shown in Figure 3.2, and initial and final deposit dimensions are given in Tables 
3.1 and 3.2. 


slide 

number 

V 

(km 3 ) 

horizontal 
coverage (km 2 ) 

L r 

(km) 

L f 

(km) 

W 

(km) 

a, 

o 

1 

0.283 

15 

2.5 

9.8 

1.7 

15.3 

2 

2100 

1500 

5.5 

79 

29 

25.9 

3 

1700 

1600 

5.5 

72 

26 

24.4 

4 

67 

1300 

3.8 

51 

26 

32.5 

5 

14 

170 

6.7 

26 

7.1 

26.4 

6 

460 

2000 

13 

58 

29 

16.4 

7 

1300 

1000 

9.0 

60 

33 

22.0 

8 

36,000 

420,000 

160 

690 

500 

25.9 

9 

0.200 

21 

5.0 

8.8 

4.2 

24.8 


Table 3.2 Landslide deposit dimensions. Values for the terrestrial slide (1) were 
calculated from maps and cross-sections from Shreve (1968), values for martian 
slides were measured from MOLA data, and values for the lunar slide (9) were 
calculated from data and cross-sections from Howard (1973). 
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Figure 3.2(b) Landslides 6 through 9. Images of Valles Marineris slides (in both parts 
of this figure) are from the Viking missions, while the Blackhawk image (Slide 1) is 
from the U.S. Department of Agriculture, the Olympus Mons aureole image (Slide 9) 
is a bitmap generated from 1/32 by 1/32 degree MOLA data (brightness is 
proportional to elevation), and the Apollo 17 landslide image (Slide 9) is from 
NASA. 
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Slide 1 . The California Blackhawk landslide (Figure 3.2a) is a terrestrial example 
of a long-runout landslide (Shreve 1968). With V= 0.28 km 3 , it is the smallest slide 
modeled here. Detailed profiles (Shreve 1968) provided useful information on the 
shape of the runout path. The lack of an abrupt change in angle between canyon wall 
and floor typical of Valles Marineris slides makes a wedge-shaped initial mass more 
appropriate for the Blackhawk slide than the block shape of Figure 3. 1 . 

The general fluidization and Bingham models produce equally satisfactory results 
(Figure 3.3). The rock mass does not fail frictionally when dry and thus the starting 
conditions for acoustic fluidization are not met. The factor of safety (FS), defined as 
the ratio of resisting to driving forces, has an average value of 1.33 (failure is 
expected to occur at FS < 1 ), although local values of FS may be lower (Zaruba and 
Mencl 1982). Therefore the slide may be considered only marginally stable for a 
friction angle of 20°. A more appropriate measure of failure probability in the context 
of acoustic fluidization is the highest friction angle that still produces movement at 
the toe of the initial shape: 1 8.9° for this model, probably close enough to the baseline 
value (given its uncertainty) that acoustic fluidization cannot be ruled out altogether if 
a suitable explanation can be found for such a low friction angle. 
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Figure 3.3 DAN results for Slide 1, the Blackhawk landslide. In each panel, the 
dashed line is the final profile (measured from Shreve ( 1 968) data) and the thin black 
line is the inferred runout path. The heavy line represents the DAN results using the 
rheology indicated. The short line at right angles to the runout path in the general 
fluidization panel marks the position at which the frictional rheology switched over to 
power law. 

A pore-pressure coefficient of 0.27 corresponds to hydrostatic pressure (with rock 
density of 2700 kg.m' 3 and water density of 1000 kg.m°), so the values of 0.65 and 
0.31 (Table 3.3) required for the frictional and general-fluidization models, 
respectively, are superhydrostatic. Experiments conducted by Iverson etal. (1997) 
show that on failure, pore pressures may increase by as much as 100 % due to soil 
compaction. Such an increase raises the pore-pressure coefficient to 0.42, above that 
required for the general-fluidization model that best matches this slide. Fluid at 
lithostatic pressure has r u = 0.5, so values in excess of this limit are physically 
impossible, ruling out the viability of purely frictional dynamics for this slide. 
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slide 

Frictional 

Bingham 


General 

Fluidization 


Acoustic 

Fluidization 

number 

r u 

To 

(kPa) 

F 

(kPa.s) 

To 

(kPa) 

F 

(kPa.s) 

r u 

To 

(kPa) 

F 

(kPa.s) 

1 

0.65 

24.5 

10 

13.5 

5.0 

0.31 

- 

- 

2 

0.56 

870 

200 

135 

20 

0.49 

245 

20 

3 

0.47 

790 

50 

49.5 

10 

0.35 

157 

20 

4 

0.58 

265 

40 

79.0 

50 

0.29 

- 

- 

5 

0.58 

75.0 

20 

21.8 

5.0 

0.00((f) = 
22.8°) 

21.8* 

5.0* 

6 

0.71 

170 

40 

65.5 

10 

0.48 

- 

- 

7 

0.58 

525 

200 

145 

50 

0.41 

146 

50 

8 

0.94 

50.0 

20 

49.2 

5.0 

0.98 

- 

- 

9 

0.41 

0.50 

0.30 

1.8 

0.50 

0.00 

1.8* 

0.50* 


Table 3,3 Summary of DAN results. Rheologies other than general fluidization 
produced unequivocal best-fits for only 2 of the 9 slides, and these are indicated by 
asterisks. The frictional phases of all slides had $ - 20°, except where indicated, and 
for acoustic fluidization r u = 0 in all cases. Where no data for acoustic fluidization are 
given, this rheology did not cause the initial mass to fail. 


we define rear runout L r as the minimum x-coordinate of the final deposit 
configuration and front runout Lf as the maximum x-coordinate. The rear runout of all 
Bingham slides is zero because the tapered upper end of the initial mass is, in general, 
too thin to produce stresses above the yield value. The relative motion of the lower 
portion of the initial mass thus produces a thin, stretched tail. In DAN, this causes the 
rear elements of the slide to lengthen significantly and the tail may not always be 
smoothly resolved - models are run with higher resolutions where necessary. 


Slide 2. One of the largest landslides considered here, Slide 2 (Figure 3.2a) is 
situated against the north wall of Ophir Chasma. Three distinct deposits are visible 
below this wall and Slide 2 is the furthest west. All three deposits have clearly 
defined tails, but further along their lengths they overlap and are difficult to 
distinguish. They are, nonetheless, good examples of unconfined long runout 
landslides. Slide 2 has a front runout of 79 km, a distance ten times that of the initial 
height of the failed material (Table 3.2). 
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The frictional and general fluidization rheologies provide comparable fits to the 
final configuration (Figure 3.4). General fluidization produces better toe and tail 
shapes, but is less accurate between x = 20 and 40 km. Required pore pressures in 
both models are high (r„ = 0.49 and 0.56 for the frictional and general-fluidization, 
respectively); the small difference between the two does, however, span the lithostatic 
value, perhaps marginally favoring the general fluidization model. Note that the 
position of the transition from frictional to power law for general fluidization was 
chosen in this case to coincide with an abrupt decrease in profile slope. The transition 
position for the acoustic fluidization rheology is, however, further back along the 
runout path since the initial frictional stage is less energetic without pore pressure and 
reaches its maximum velocity sooner. 
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Figure 3.4 DAN results for Slide 2. The observed final profile (dashed line) is 
measured from MOLA data. The dotted line is the initial shape of the landslide mass. 
Note that the runout path and initial shape have rounded comers, a result of DAN’s 
spline routine which generates only smooth curves. Acoustic fluidization is a special 
case of general fluidization (pore pressure is zero in the frictional stage). 


Slide 3. This landslide is the furthest east of the North Ophir Chasma slides 
(Figure 3.2a). Like Slide 2, the scar in the wall of the chasm is clearly contrasted with 
its surroundings, as is the tail of the deposit, and most deposit dimensions are similar 
to those of Slide 2 (see Table 3.2). Although both the frictional and general- 
fluidization rheologies match much of the profile (Figure 3.5; x = 20 - 55 km), only 
the latter matches the slide tail and the long, thin runout to* = 72 km. Furthermore, 
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significantly lower pore pressure is required for general fluidization compared to 
friction ( r u = 0.35 vs. 0.47). Note that the position of the transition from frictional to 
power law in the general-fluidization model was also chosen here to coincide with an 
abrupt decrease in profile slope. 



Figure 3.5 DAN results for Slide 3. 


Slide 4. This relatively small landslide (number 1 in Lucchitta 1979) is located 
against the south wall of Gangis Chasma. It is the only slide of those considered here 
for which MGS Mars Orbiter Camera (MOC) data are available. An image covering 
the length of the slide (Figure 3.6) reveals detailed structure, including the scar 



features described above. The deposit has longitudinal ridges and grooves, and 
transverse ripple-like structures reminiscent of a turbulent fluid. 
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Figure 3.6 Two parts of a MOC image (M07/03026) oriented along the length of 
Slide 4. Each part is just under a tenth the length of the original image (whose 
footprint is visible in Figure 2a). Image width covers approximately 1 .5 km. (a) 
Landslide reentrant with dark upper layer and light talus slopes below, (b) Transverse, 
ripple like structures reminiscent of turbulent flow. Although the contrast in the entire 
image has been stretched, the two strips have not been stretched relative to each other. 


General fluidization produces the best fit (Figure 3.7) with only modest pore 
pressure (Table 3.1). The Bingham fluid reproduces the bulbous toe shape, but fails to 
mimic the tail. Because the structure of the runout path is somewhat simplified, the 
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apparent deepening of the deposit near the toe might also be interpreted as a constant 
thickness deposit overlying a bulge in the runout path. Despite the steep initial slope 
of the runout path, acoustic fluidization does not initiate failure. This is because much 
of the material lies initially over the level section of the runout path and its weight 
does not contribute to the downslope driving force. The highest friction angle that 
produces motion in the toe of the initial shape is 14.5°, significantly less than the 
baseline value. 



x [km] 


Figure 3.7 DAN results for Slide 4. 

Since the initial and final cross-sectional areas are very similar (Table 3.1), the 
decrease in area caused by the great degree of widening of the slide (Figure 3.2a) has 
to be compensated by a high (10 7 m 3 .m'') entrainment rate, which may imply that the 
canyon wall in this region is weaker than that of the other slides. 
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Slide 5. Situated against the north wall of the eastern reaches of Gangis Chasma, 
Slide 5 (number 8 in Lucchitta 1979) is the smallest Martian landslide modeled (V = 
14 km 3 ). The shape of the deposit toe is influenced by its contact with another slide 
emanating from the opposite canyon wall (Figure 3.2b and 3.8). Consequently it is 
not surprising that none of the rheologies, except perhaps the Bingham fluid (which 
was run with 30 slide elements to improve tail resolution), reproduces the narrow 
hump at the deposit toe. Overall, the acoustic fluidization rheology provides a 
comparably accurate fit; note that the friction angle had to be increased to 22.8° from 
the baseline value of 20°. The initial height H (Table 3.1) of the failed material is not 
significantly less than that of other Martian landslides, but because the deposit 
volume is so small, the depth D of the deposit (Figure 3.1; Table 3.1) is also small. 
This gives rise to a greater rear runout, hence the reduced need for pore pressure that 
allows acoustic fluidization to model this slide successfully. In the available Viking 
images of Slide 5 no obvious scar is visible on the canyon wall and H was thus 
chosen to correspond to the full height of the wall, a rule of thumb supported by 
observations of most other large scale Valles Marineris slides. 
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Figure 3.8 DAN results for Slide 5. The best-fit was acoustic fluidization: 
introduction of pore pressure made the results less accurate, hence the coinicidence of 
optimized acoustic and general fluidization results. 
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Slide 6. This landslide (Figure 3.2b; number 9 in Lucchitta 1979) is located at the 
west end of Coprates Chasma. It is a relatively small slide, with a distinctively level, 
square tail (Figure 3.9). This may be due to the wholesale slipping of the topmost 
section of the original mass. DAN does not account for the motion of large rigid 
blocks so the main priority of Slide 6 models was to match front runout. None of the 
models produces a satisfactory fit. Acoustic fluidization does not initiate failure; the 
highest friction angle that produces motion in the toe of the initial shape is 16.3° (FS 


= 1.24). 



Figure 3.9 DAN results for Slide 6. 


Slide 7. Located against a west facing wall in the region between Ophir and 
Melas Chasma, this slide (Figure 3.2b; number 1 1 in Lucchitta 1979) is similar in 
volume and dimension to the north Ophir Chasma slides (Table 3.2), but has a shorter 
runout. General fluidization, which produces the best fit (Figure 3.10), thus requires a 
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high yield strength (145 kPa). The required pore-pressure coefficient (0.41) is 
between hydrostatic and lithostatic values. The frictional rheology also yields a good 
fit to this profile but requires superlithostatic fluid pressure ( r u = 0.56). 



Figure 3.10 DAN results for Slide 7. 

Slide 8. we ran a model of the westward facing Olympus Mons aureole deposit. 

In a bitmap generated from MOLA data (Figure 3.2b), where each picture element 
has a brightness proportional to its elevation, the deposit can be seen overlying a 
darker (and therefore lower elevation) deposit. The scarp near the tail of the deposit is 
assumed to be the landslide scar, and the initial mass is wedge-shaped in order to 
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simulate the pre-slide shield topography. Post-emplacement plains and surface 
deposits covering the rear of the aureole make it difficult to infer the deposit 
configuration in this region, and we thus assume a uniform deposit thickness. Three 
different profiles were constructed in order to get an overall picture of the aureole 
topography. The runout path dips towards Olympus Mars proximally, and away from 
it distally, agreeing with observations by Francis and Wadge (1983). Neither 
frictional, Bingham, nor general-fluidization models produce particularly good fits 
(Figure 3.1 1 ) and most notably the required pore-pressure coefficients must exceed 
0.9 to achieve the 690 km front runout. Failure is not initiated for the dry frictional 
precursor to acoustic fluidization: FS = 2. 1 0 and the highest friction angle that 
produces motion in the toe of the initial shape is 3.8°. we conclude that the Olympus 
Mons aureole is not a landslide deposit but the result of an alternative process such as 
gravity sppreading (Tanaka, 1985). 



Figure 3.11 DAN results for Slide 8, the northern Olympus Mons aureole deposit. 
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Slide 9. we modeled a large lunar landslide near the Apollo 17 landing site whose 
volume is similar to that of the Blackhawk slide (Table 3.2). Although detailed 
topographic data is unavailable, approximate rear and front runouts (5.0 and 8.8 km) 
and estimates of average proximal and distal deposit thicknesses (20 and < 10 m) are 
known (Howard 1973). Since no obvious scar is visible, we assumed that the 
landslide material initially covered the full height of the mountain massif down which 
it slid, we did not include width or volume variations in this model. The initial shape 
is very shallow: with D = 22 m the depth normal to the runout path is only 22 sin 
24.8° = 8.9 m (Table 3.1). Since the lunar regolith is 5 to 10 m deep (Quaide and 
Oberbeck 1973) it is unlikely that bedrock failure was involved. 



Figure 3.12 DAN results for Slide 9, the Apollo 17 landslide. No measured final 
profile is available, but the known required rear and front runouts are marked with 
vertical dashed lines. 
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Acoustic fluidization produces the most accurate combination of rear and front 
runouts (Figure 3.12). This conforms to the clear requirement that lunar landslide 
“fluidization” must proceed in the absence of a fluid. The introduction of even small 
pore pressures produces a longer rear runout than required. The frictional rheology 
requires substantial fluid pressure to achieve the desired front runout ( r u = 0.41), 
confirming the need for a fluidizing process to describe its dynamics. 

The small depth of the initial shape requires low yield strengths to promote 
failure. A simple sliding-block model can be used to estimate yield strength 
(e.g., McEwen 1989): 

r o = Pgd sin P 

where p is density, g is gravitational acceleration, d is the deposit toe thickness, and P 
is the slope of the runout path at the deposit toe. Applying appropriate values of these 
variables, including a slope of 5° (a reasonable upper limit) and a toe depth of 10 m, 
we obtain a yield strength upper bound of 3.8 kPa, which compares well with values 
from DAN models (Table 3.3). 

As a result, the Bingham rheology produces a deposit thinner than the estimate 
from Howard (1973) of 10 to 20 m. On the other hand, the frictional rheology 
produces a short, bunched deposit with a maximum depth of 50 m. Only acoustic 
fluidization produces the desired depths. 

we repeated the simulation with H reduced to half its original value. With the 
initial shape occupying only the lower half of the mountain slope, depth D was 
doubled to conserve volume. Despite the lower potential energy of the initial mass, 
acoustic fluidization remained the best rheology. Higher yield strengths were required 
to maintain the same runout (1.20 kPa for the Bingham rheology and 1 .27 kPa for 
acoustic fluidization). 
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3.4 Discussion and Conclusions 

3.4. 1 Best overall rheology 

The general fluidization rheology is the obvious best-fit for Slides 3 and 4 and can 
be considered a comparable best-fit for Slides 1, 2, and 7. None of the selected 
rheologies produces good fits for Slides 6 and 8, which we attribute to the effects of 
large, discrete blocks and a non-landslide origin, respectively. Acoustic fluidization is 
the best-fit for Slide 9. For Slide 5 the Bingham rheology is clearly superior to 
general fluidization, but this slide is compromised by longitudinal confinement. 

The advantages of general fluidization can be described by considering the 
shortcomings of the other rheologies. The frictional rheology has a friction angle 
(20°) close to the failure angle ai in many cases (Table 3.2). This limits the reach of 
the slide, producing small rear and front displacements. Lower friction angles would 
produce greater runouts but the current choice already constitutes a physically 
reasonable minimum. Instead, large pore pressures are required to achieve the desired 
runouts and these frequently exceed the likely maximum (corresponding to r u - 0.42 
as described above). With the exception of the fluidization rheologies in the Olympus 
Mons aureole model, only the frictional rheology requires superlithostatic pore 
pressures ( r u > 0.5), emphasizing the inability of this rheology to explain long runout 
landslides (Cleary and Campbell 1993). 

The frictional profiles are typified by a linear decrease in deposit thickness in the 
positive x direction (if pore pressure is zero, the surface of the deposit comes to rest at 
the friction angle of the material). Indeed, some of the observed deposit surfaces (e.g. 
Slides 2, 3, and 7) have linear sections, suggesting that the governing rheology does, 
at times, exhibit frictional behavior, pointing towards the use of two-stage fluidization 
processes such as those considered here. 
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The Bingham rheology consistently fails to produce the correct tail shape, even 
when model resolution is increased. Toes are frequently too deep, indicating that the 
yield strength is greater than values calculated directly from deposit thickness 
(McEwen 1989). Actual yield strengths are likely to be smaller, probably similar to 
the power law yield strengths obtained by the general fluidization rheology. 

The general fluidization parameter values are mostly reasonable, although yield 
strength and apparent viscosity values are difficult to judge, being poorly constrained 
(viscosity is probably between 10 and 1000 kPa.s in dry sand, Goetz and Melosh 
1980). The pore-pressure coefficient has an average value of 0.40, a physically 
plausible value given additional pressurization due to soil compaction (Iverson et al 
1997). There are at least four potential scenarios in which geologically recent martian 
landslide formations were lubricated: 

1. As noted in the introduction, pore pressures in DAN models refer to conditions 
at the failure surface and are not representative of the entire deposit. Thus, for a fluid 
to promote greater runouts, it is not necessary for the initial material to be saturated at 
all depths, especially for power-law fluidization rheologies in which shear is 
concentrated in a basal layer. The top 300 to 500 m of near equatorial crust is 
expected to have been dessicated over the duration of current atmospheric conditions 
(Fanale et al 1986, Clifford 1993). Since most martian slides have an initial depth of 
a few kilometers (Table 3.1), a substantial part of their mass may have been saturated 
with ice at the time of failure. The frictionally dissipated energy produced during the 
slide event is, however, unlikely to be sufficient to melt ice for lubrication. The 
energy produced per unit time, per unit mass is vgs\n(a\)FS where v is the depth- 
averaged downslope velocity of the slide (Iverson et al 1997). we integrated this 
energy over time for each slide element in dry frictional DAN models and found that 
on the order of 1 kJ/kg is produced by material on the initial, steep section of the 
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runout path. For ice at 0 °C, 334 kJ/kg is required for melting, so the supply is two 
orders of magnitude too small. 

2. Alternatively, water may have been present in the liquid phase prior to failure. 
A cryosphere confining an underlying saturated aquifer, as suggested by Carr (1979) 
and Clifford (1993), might have a base shallow enough to fall above the landslide 
failure zone. Reasonable present-day values of crustal thermal conductivity (2.0 W.m' 
’.K' 1 ) and geothermal heat flux (30 mW.m' 2 ) suggest a melting point depth of 2.3 km 
at the equator, assuming that melting occurs at 252 K, a likely value given pressure 
and solute effects (Clifford, 1993). Higher heat flux estimates (e.g. 45 mW.m' 2 from 
Toksoz 1978) yield depths as low as 1.5 km, within the depth of most Valles 
Marineris landslides. Notably, of the Valles Marineris slides modeled, only Slide 5 
has a depth less than 3 km, and acoustic fluidization is its most successful rheology. 

The recent formation of a near-surface water-related gully feature at only 27° 
latitude (Malin and Edgett 2000) supports the idea of a water table only a few 
hundred meters deep, but this feature does not preclude the existence of an underlying 
cryosphere. Other sapping features have been placed contemporaneously with the 
landslides (Peulvast et al. 2001), and near-equatorial outflow channel drainage may 
have continued well into the Amazonian (Gulick 2001), supporting the idea of water 
in the Tharsis region at the time of landsliding. 

3. If lakes existed in the Valles Marineris (Carr 1996, Lucchitta et al. 1994), water 
may have infiltrated into the walls, providing sufficient pore pressure for later 
landslide fluidization. Although our principal argument here is for the presence of a 
fluid for lubrication purposes, water in this scenario might also have been partially 
responsible for the failure itself. If the Valles Marineris drained in a relatively short 
period (e.g., via the catastrophic outflow channels in the case of open depressions, 
and by flow through permeable aquifers in the case of closed depressions) leaving 
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lacustrine deposits which correspond (at least in part) to the interior layered deposits 
of the present day, then the newly exposed canyon walls, no longer under pressure 
from a large body of water and still highly saturated, would be unstable, resulting in 
widespread landsliding. Lucchitta et al. (1994) suggested that the Valles Marineris 
spur-and-gully formations may have formed subaqueously. They would thus have 
formed prior to the landslides while the lakes were still in place. This scenario is then 
in agreement with the observed relative ages of the landslides, spurs and gullies, and 
the interior layered deposits. Problems with this means of landslide emplacement 
include the ability of the crust to hold water for a sufficiently long period. With a 
length scale on the order of 10 km, and a permeability of 10‘ 17 m 2 (a likely crustal 
lower bound, Harrison and Grimm 2002), a 1 myr period would be available for 
water-lubricated landslide formation (this is the time scale for diffusive drainage of 
fluid from an initially saturated aquifer 10 km thick). This is a reasonable lower 
bound only, and it is likely that landsliding occurred over a longer period (given the 
sometimes large differences in deposit age, Lucchitta 1979), requiring unreasonably 
low permeabilities. 

4. In the absence of water, CO 2 may have been a suitable fluidizing agent. 
Hoffman (2000) proposed a “white Mars” in which CO 2 is the principal fluid at work 
in the catastrophic flood channels which open into the northern highlands. 
Depressurization of subsurface CO 2 liquid causes it to flash into the gaseous phase 
and mix broken landslide material into a turbulent cloud of dust, rocks, ice, and gas 
(Hoffman, 2000). Continual degassing of liquid from the landslide mass provides an 
ongoing source of lubrication. Hsu (1975) simulated the terrestrial Elm landslide with 
a mixture of silt and dry ice whose degassing did indeed produce long runouts. 

we conclude this section by mentioning some aspects of landslide dynamics that 
may contribute to long runouts but are not simulated in DAN. The rear of many of the 
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Martian deposits contain unbroken segments of the original landslide mass that have 
slumped or rolled down the runout path, possibly producing greater rear runout than 
material composed of uniformly sized small-grained particles. In most cases, these 
blocks are not sufficiently isolated from smaller grained particles to justify their 
exclusion from DAN models, with the exception of Slide 6 whose tail has a large 
level section, indicating that the upper part of the initial mass slipped in one piece 
down the failure plane. The dynamical influence of large blocks, and possibly other 
phenomena, may mean that actual pore pressures are somewhat less than those 
predicted by DAN models, although they are likely to remain non-zero. This effect 
will be offset by the increase in required pore pressures if higher friction angles are 
found to be more appropriate. 


3.4.2 Other Rheologies 

The Bingham rheology produces the best results for Slide 1. A yield strength of 
24.5 kPa and a viscosity of 10 kPa.s are, however, significantly smaller than the 
values used to model rock avalanches (300 kPa and 40 kPa.s respectively, Hungr 
1995). 

Acoustic fluidization produces the best results for slides 5 and 9. Although these 
slides are unique in having relatively small values of D, models of Slide 9 with twice 
the value of D (and half the value of H) yield the same rheology ranking, indicating 
that the suitability of acoustic fluidization is determined by other factors such as 
runout slope. 

A geometrical feature of significance in Slides 4 and 8 (and to a lesser extent in 
slides 2, 3, and 7) is a large initial “footprint”. The part of the initial mass which rests 
on the horizontal canyon floor does not contribute to the downhill force that promotes 
movement, leading to shorter runouts which must be compensated by higher pore 
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pressures or lower yield strengths. Although Slide 4 has a large footprint, its runout 
path is steeper than that of the other slides, resulting in a relatively good performance 
by acoustic fluidization. 

The long runout of the Olympus Mons aureole is comparable only to submarine 
landslides, such as those near Hawaii (Moore 1964) which depend on a high degree 
of lubrication. The shoreline of a northern hemisphere ocean, if one existed at the 
time of failure, would probably have had a mean elevation of about -3800 m (Head at 
al 1999) with excursions as high as -2000 m at Tharsis (due to subsequent uplift of 
the region), which only just coincides with the lowest parts of the aureole deposit. 
Complete submersion would have required sea levels as high as 7000 m. The 
improbability of submarine formation and the extreme values of many of the DAN 
parameters suggest that the mode of emplacement was something other than 
catastrophic landsliding. Favorable alternatives include ice-lubricated gravity 
spreading (Tanaka 1985), a low strain-rate process that may have acted over a period 
as long as 1 myr. Spreading may have occurred predominantly along a decoupling 
zone lubricated by ice that possibly originated from earlier lateral spreading of 
hydrothermal waters driven by volcanic activity within Olympus Mons (Tanaka 
1985). 


3.4.3 Gravity Effects 

The ratio HI Lf 4s a useful, although not perfect, measure of the extent to which a 
landslide runout is “long”. McEwen (1989) and Hsu (1975) provide HI L/\s. volume 
data for Mars and the earth, respectively (Figure 3.13). H I Lf decreases with volume, 
meaning that larger landslides have longer runouts relative to their initial height. Also 
of note is the offset between the Martian and terrestrial trends, evidence that the 
rheology governing the slides cannot be purely frictional. For a given set of parameter 
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values, frictional slides produce a constant HI Lf value independent of V and 
gravitational acceleration because the basal resistance scales with gravity and volume 
and therefore with the downslope acceleration. For the Bingham rheology and the 
power law part of the fluidization rheologies, the basal resistance does not scale with 
the downslope acceleration because it is based on a gravity- and volume-independent 
yield strength. The result is a dependence of/// Lf on volume and the emergence of 
separate HI Lf trends for different gravity environments. 



Figure 3.13 Observed trends in ratio HI Z/vs. V for Valles Marineris {McEwen 
(1989)) and terrestrial and lunar landslides ( Hsu (1975)). Values for Slides 1 to 9 are 
shown with numbered circles. Slides 2 and 3 correspond to two of McEwen’s data 
points but do not coincide exactly with them because of differing volume 
measurements. Curves generated by simplified DAN models of Slide 5 with varying 
volume under both terrestrial and martian gravity demonstrate the ability of the 
Bingham (B) and general fluidization (G) rheologies to mimic the observed trends. 


McEwen (1989) recognized the need for a gravity independent yield strength (see 
also Dade and Huppert 1998) and in Figure 3.13 his point is demonstrated by a set of 
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DAN models of varying volume, run with both the Bingham and general fluidization 
rheologies in martian and terrestrial gravity. A uniform width, zero entrainment 
version of the Slide 5 DAN model was scaled spatially to obtain slides of different 
volumes but identical aspect ratio ( D/H ). The rheological parameter values that 
produced the required front runout in the original scale were used for all other scales 
on both planets, reflecting the assumption that yield strength and viscosity are 
independent of gravitational acceleration. Slide 5 has a fairly low yield strength: 
experiments with higher yield strength slides (such as Slide 2) show that the curves in 
Figure 3.12 shift upwards together without significantly changing their offset. Note 
that the symbol representing Slide 5 (Figure 3.13) falls over the intersection of the 
two solid curves at V = 14 km 3 , implying that both sets of rheological parameters 
were chosen so that Lf, and therefore H / Lf, coincide at the original scale. 

The Bingham and fluidization rheologies reproduce the first-order features of the 
observed trends extremely well. The general fluidization curves level out with 
decreasing volume because the efficiency of the power law stage decreases (the yield 
strength remains constant while the thickness of the deposit decreases), and the 
frictional stage becomes the dominant influence on runout. This departure reflects the 
fact that fluidization is unlikely to occur at low volumes (relative to a specific 
gravitational acceleration). The suitability of the Bingham rheology is also somewhat 
restricted: it does not initiate failure for volumes below about 0.2 km 3 for Martian 
slides and 0.001 km 3 for terrestrial slides. It must be remembered that all the curves 
in Figure 3.13 were generated with one set of rheological parameter values. It is 
physically unlikely, however, that every martian landslide can be described by the 
same parameter values. There will be a dependence on grain size, rock or soil density, 
fracture distribution, and other local phenomena. Parameter variations are also 
evident in the DAN results, although these may be exacerbated by the simplified 
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nature of the model that excludes small scale, unobservable variations in failure 
geometry, and other second order effects. Overall, Figure 3.13 provides further 
support for McEwen’s claim of a gravity independent yield strength, and extends it to 
general and acoustic fluidization, processes which include the frictional rheology 
(which alone cannot account for the observed trends). 

The west Olympus Mons aureole (Slide 8, Figure 3.13) lies well below the Bingham 
and general fluidization trends for Mars. This reflects its unusually long runout 
(hence low HI Lj) which our rheologies were only able to reproduce with high pore 
pressures and/or low yield strengths. The Blackhawk landslide (Slide 1) is noticeable 
as the only data point lying close to the terrestrial trend, as expected. Its position 
relative to the curves generated by Slide 5 parameter values is a combination of lower 
yield strength (which lowers HI Lj) and higher aspect ratio (which raises HI Lj) than 
Slide 5. The Apollo 17 landslide (Slide 9) lies between the terrestrial and martian 
trends in Figure 3.13. Because the moon has the lowest gravitational acceleration of 
the three bodies, its landslides might be expected to lie above the martian trend 
(vertically separated by about 2.3, the ratio of martian and lunar gravities); however 
Slide 9 falls below this trend. At least two factors may have affected its value of H / 
Lf\ 1) impacts responsible for the initial failure may have contributed energy to the 
subsequent motion, resulting in longer runouts than expected (this may have affected 
other lunar slides close to impact craters, such as the Tsiolkovsky landslide (Wu etal. 
1972, Guest 1971, Figure 3.13), and 2) the slide is probably an avalanche of regolith 
material (which is as deep as 15 m in places, Quaide and Oberbeck 1973) and not a 
bedrock landslide. 
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CHAPTER 4 
OUTFLOW CHANNELS 

4.1. Introduction 

On Mars, the best evidence for groundwater flow probably lies in the outflow 
channels, structures larger than any terrestrial flood-related features, with lengths and 
widths on the order of hundreds of kilometers (Sharp and Malin, 1975; Carr, 1996) 
Two-thirds of all outflow channels are located around the Tharsis rise, with almost all 
of these debauching into Chryse Planitia (19/30 and 16/19, respectively, using the 
tabulation of Clifford and Parker, 2001). The source-region structure of the channels 
suggests the release of groundwater following the local disruption of a frozen 
cryosphere (Carr, 1 979). Circum-Chryse outflow formation began some time after the 
end of heavy bombardment 3.8 Gyr ago (Tanaka, 1992; Carr, 1996) continuing to the 
end of the Hesperian 2.9 to 3.3 Gyr ago (Hartmann and Neukum, 2001; see Clifford 
and Parker, 2001). Flow may have been catastrophic (Baker and Milton, 1974), with 
flooding lasting from a few days (e.g. Carr, 1996) to years (Williams and Phillips, 
2000; de Hon et al., 2003), possibly in several episodes (Scott and Dohm, 1990; 
Williams and Phillips, 2000). 

A fully consistent description of outflow channel formation remains elusive 
primarily because of the difficulty in explaining the large volumes of water involved 
and the required flow rates. At least 6 x 10 6 km J (Carr, 1996) are thought to have 
passed through the circum-Chryse channels; however, such estimates require 
specification of the volume of water required to remove a unit volume of sediment 
(the sediment- to- water ratio), but this parameter is poorly constrained. Komar (1980) 
suggests values as high as 0.67 for the large outflows such as Kasei, whereas de Hon 
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et al. (2003) suggest values as low as 10" 4 for smaller channels such as Maumee and 
Vedra Valles. 

Estimated flow rates range from 3 x 10 6 m 3 /s (Williams et al., 2000) to 10 10 m 3 /s 
(Carr, 1996; Mosangini and Komatsu, 1998). The lower values in this range might be 
explained by Darcy flow through a porous medium under high pressures. Carr (1979) 
considered the disruption of a cryosphere due to lithostatic pore pressures, showing 
that permeabilities as high as 3 x 10 9 m 2 are required to achieve flow rates of 10 6 m.s' 
'. Such permeabilities are comparable to the highest measured terrestrial values 
(Robertson, 1 974) and are about an order of magnitude greater than typical high 
permeability rocks such as Hawaiian volcanics (Davis, 1969). Because Darcy flow 
through a porous medium of moderate permeability and porosity might not be able to 
sustain such high flow rates (Head et al., 2003), water drawn from distal parts of the 
aquifer might have to be stored in a local, highly permeable “buffer” or reservoir 
ready to deliver high flow rates. The most likely candidates are the Valles Marineris, 
which may have held ice-covered lakes (Squyres, 1989; Lucchitta et al., 1994), 
possibly containing enough water to explain the required discharge volumes. 
Alternatively, the detailed analysis of channel geometry with high resolution images 
(e.g. Williams and Phillips, 2003) may continue to suggest lower flow rates than 
previously thought, obviating the need for a local reservoir. 

While the local geologic and hydrologic properties of the outflow source regions 
are complex and still largely unknown, the regional-to-global aspects of outflow 
formation - which are controlled by the interaction between groundwater flow, 
storage, recharge, and discharge - are easier to address. In particular, we look at the 
role of groundwater flow in circum-Chryse outflow channel formation, starting with 
the work of Clifford (1993), and Clifford and Parker (2001), henceforth named C&P, 
who present a hydrologic model explaining outflow channel and south polar cap 



77 


formation. The model starts during the Noachian with a relatively warm atmosphere 
and high geothermal heat flux (95 mW.nT 2 , C&P), conditions that allow unconfined 
groundwater flow. The inevitable result, according to C&P, is net groundwater flow 
toward the low elevations of the northern hemisphere, where an ocean develops. 
During the late Noachian the atmosphere begins to thin and cool, the northern ocean 
starts to freeze, and ice sublimates from its surface, moving to the south pole where it 
condenses as ice onto a polar cap. The subsequent growth of the cap causes the local 
melting isotherm to rise until it intersects the crustal surface. At this stage, the base 
of the cap begins to melt and water seeps into the crust, forming a groundwater 
mound which drives flow toward lower latitudes, raising global water table levels and 
providing a source for the outflow channels. Meanwhile, the cooling atmosphere is 
accompanied by the development of a cryosphere which becomes thick enough to 
confine the global aquifer and, when locally disrupted, gives rise to catastrophic 
outflows. 

Here we defer analysis of Noachian conditions and focus on Hesperian 
groundwater dynamics. We test the C&P south pole recharge model and an 
alternative source of groundwater recharge and hydraulic head for the circum-Chryse 
outflows: Tharsis. 

Tharsis has been the primary center of volcanism throughout Martian history and 
many workers argue for a causal relationship between Tharsis volcanism and outflow 
channel formation (Masursky, 1986; Tanaka and Chapman, 1990; Head and Wilson, 
2001; Chapman and Tanaka, 2002). Heat may have thinned the cryosphere (Tanaka et 
al ., 2002), perhaps melting it through in places (Wilson and Head, 2002). McKenzie 
and Nimmo (1998) go further to suggest that dike intrusions were responsible for the 
formation of the Valles Marineris and attendant cryosphere melting and outflow 
discharge. 
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The role of equatorial volatiles in outflow formation extends to the atmosphere 
and centers around variations in the planetary obliquity (tilt of rotation axis relative to 
the ecliptic). The current value of 25° leads to average annual temperatures about 60 
K warmer in the tropics than at the poles. However, the obliquity varies periodically 
on a scale of about 0. 1 Myr (Ward, 1 992) and is thought to have reached values as 
high as 60° (Laskar and Robutel, 1993; Touma and Wisdom, 1993). The effects of 
high obliquity on climate are profound and include the development of a cold 
equatorial zone in which surface ice is perennially stable (Jakosky and Carr, 1 985; 
Mischna et al., 2003; Haberle et al., 2003; Richardson and Wilson, 2002) . It has 
been suggested that equatorial ice could, at high obliquity, be sourced by water 
sublimated from the poles (Jakosky and Carr, 1985; Coleman and Dinwiddie, 2003; 
Coleman et al., 2003; Mischna etal., 2003; Mischna, personal communication; 

Mellon et al., 1997; Haberle et al., 2000) and be preferentially deposited over regions 
of high topographic elevation (Mischna et al., 2003) in analogy to the earth. The low 
latitude and high elevation of the Tharsis rise make it a good candidate for 
precipitation at high obliquity. Indeed, much evidence of mid- to high-latitude ice 
remains today, including shallow subsurface traces detected by the Mars Odyssey 
Gamma Ray Spectrometer (Boynton et al., 2002), ice-rich layered deposits from 
recent ice ages (Mustard etal., 2001 and 2003), and gullies possibly formed by 
melting snow transported from polar to mid-latitudes during high obliquity 
(Christensen, 2003). Ice-related features are also observed at low-to mid-latitudes 
latitudes and include rootless cones indicative of volcano-ice interactions (Greeley 
and Fagents, 2001), rampart craters (Allen, 1979), and glacial landforms flanking 
volcanic edifices (Head and Marchant, 2003). 

In the sections ahead we bring to light problematic aspects of the C&P model, and 
we infer from the evidence given above that groundwater recharge over Tharsis 
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constitutes a preferable alternative to South Pole recharge. We support our argument 
with numerical simulations that test and compare both models. 

4.2. Model 

4.2.1 MODFLOW-2000 

We simulate saturated groundwater flow with MODFLOW-2000 (Harbaugh et 
al., 2000), a recent version of the well-established model of the United States 
Geological Survey. The code performs three-dimensional transient groundwater 
simulations, solving the governing equations for porous media flow on a discrete 
finite-difference grid. Because thermal effects are ignored, a complete description of 
the behavior of the system is possible with a single partial differential equation in the 
hydraulic head h: 


V.(K .Vh)+Q = S — (4.1) 

s dt 

where K is the hydraulic conductivity tensor, Q comprises sources and sinks, S s is the 
specific storage of the porous medium, and t is time. The hydraulic conductivity, 
which in the models of this work is generally heterogeneous (but isotropic), is related 
to other properties of the porous medium and its fluid as follows 

K=-^ (4.2) 

M 


where k is the permeability of the medium, p is the fluid density (1000 kg.m' J ), g is 
the Martian gravitational acceleration (3.71 m.s' 2 ), and p is the dynamic viscosity of 
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the fluid (10‘ 3 Pa.s). The source term Q consists of surface recharge and drainage, 
described in detail below. 

Permeability, while isotropic, varies with depth according to globally assimilated 
terrestrial measurements (Manning and Ingebritsen, 1999), and is scaled to the 
Martian environment by multiplying depth by the ratio of Martian and terrestrial 
gravitational accelerations (C&P, g^/ gE = 3.7/9.8): 

log* = - 12.65 -3.21ogz (4.3) 


where z is the depth in km. 

We consider two porosity relationships. The first, used by C&P, is Athy’s Law 
(1930) scaled from lunar seismic data (Binder and Lange, 1980): 

n - n 0 exp(-z/2.82) . (4.4) 

where no is the surface porosity. 

The second relationship is based on the Kozeny-Carman equation (Ingebritsen and 
Sanford, 1998) which relates porosity and permeability. We use the following 
version modified for vesicular basalts (Saar and Manga, 1999): 

k = An 4 

where A = 8.9 x 10' 7 m 2 contains information regarding the tortuosity, shape, and 
specific surface area of vesicular basalt porosity. Using the porosity-depth 
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relationship of Equation 4.3, we have 

log n = - 1 .650 - 0.8 log z . (4.5) 

The specific storage coefficient S s in Equation 4.1 controls, under confined 
conditions, the volume of water released from the aquifer due to compression of the 
rock matrix per unit increase in hydraulic head. For water released under unconfined 
conditions, the porosity of the rock must be taken into account and the storage 
formulation is best expressed in terms of the storativity: 

S = S y +bS,. (4.6) 

Here b is the thickness of the saturated portion of the aquifer and S v is the specific 
yield (that part of the pore space that can be drained, and therefore approximately 
equal to the porosity). Commonly the specific yield is used on its own under 
unconfined conditions (Fetter, 1994) since it is generally much larger than bS s - this is 
the approach MODFLOW-2000 adopts. 

4.2.2 Boundary Conditions 

The elevation of the upper model boundary corresponds to the roof of the Martian 
aquifer, and depends both on the topography and the cryosphere thickness. The 
topography is taken from Mars Global Surveyor MOLA data and therefore assumed 
to be that of Mars at present, with minor modifications. This assumption can be cast 
as a statistical argument, that the hypsogram of Mars during the past is approximately 
that of the present (C&P). However, there is direct evidence that the bulk of the 
Tharsis rise was in place by the end of the Noachian (Phillips et ah, 2001), and 
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therefore the overall topography of Mars has changed very little since then. We 
remove the large Amazonian shield volcanoes (Olympus Mons, Tharsis Montes, Alba 
Patera, Elysium Mons) because their high elevations would artificially influence 
Hesperian hydraulic heads. Crvosphere thickness is determined by the surface 
temperatures and geothermal heat flux being simulated, and on an assumed melting 
temperature of 252 K (C&P). The lower model boundary is 8.5 km below the 
surface, corresponding to the depth of pore space compaction assumed by C&P, 
determined by a porosity of 1% in Equation 4.4. The alternative relationship of 
Equation 4.5, however, falls more steeply than Equation 4.4, reaching 1% at 2.7 km, 
so for models using this porosity relationship we relax the somewhat conservative 
threshold of C&P, maintaining a depth of 8.5 km, corresponding to a compaction 
porosity of n — 0.4%. Integration of hydraulic conductivity with depth in an aquifer 
with a 3 km thick cryosphere (the average thickness for latitudes between -30° and 
30° in model 2) shows that 91% of the flow in the tropics is expected to occur in the 
top 8.5 km of crust, so the flow here will not be affected by the presence of the 
compaction threshold. At the poles, the thicker cryosphere results in a smaller fraction 
(50%) of the total flow above a depth of 8.5 km, and the compaction threshold will 
thus have a moderate effect on the local flow. 

All models have six layers whose thicknesses increase with depth, allowing the 
steep gradient of the permeability and porosity at shallow depths to be modeled 
accurately while keeping the number of layers, and thus the computational overhead, 
small. The models are cast in an equal-area sinusoidal projection to conserve mass 
and reduce distortion at the poles. Areal coverage of each finite difference cell is 
2800 km 2 . 

The left and right model boundaries follow a longitude of 180° (East) along which 
the only major topographical variation is the dichotomy boundary whose steepest 
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gradient is oriented approximately S-N, making this longitude a suitable artificial no- 
flow boundary as flow lines will tend to align with it. 

We assume that the Valles Marineris crust is in thermal equilibrium with the 
atmosphere and therefore has a cryosphere thickness similar to that of the nearby 
crust, i.e. we do not allow leakage of water through the Valles Marineris walls to 
produce lakes. The model is transient and run for 900 Myr, corresponding to the 
period between the end of heavy bombardment (3.8 Gyr, Tanaka et al., 1992) and the 
lower bound on the Hesperian-Amazonian boundary (2.9 Gyr, Hartmann and 


Neukum, 2001). 

Initial head 

is taken to 

correspond 

to the 

aquifer surface, 

representing an abrupt end to 

the Noachian 

with instantaneous 

freezing of the 

cryosphere. 






Source Region 

Source 
Volume 
10 4 km 3 

Associated 

Valleys 

Additional 
Volume 
10 4 km 3 

Total 
Volume 
10 4 km 3 

Total Water 
Volume* 
10 4 km 3 

I. Aram Chaos 

2.8 

- 

- 

2.8 

4.20 

2. Aromatum 
Chaos 

4.5 

Shabaltana 

- 

4.5 

6.75 

3. Aureum- 
Arsinoes Chaos 

13.2 

Simud, Tiu 

27.2 

40.4 

60.6 

4. Echus Chasma 

53.5 

Kasei 

36.9 

90.4 

136 

5. Hydaspis Chaos 

4.5 

Tiu 

8.7 

13.2 

19.8 

6. Hydraotes 
Chaos 

11.7 

Simud, Tiu 

11.7 

23.4 

35.1 

7. lani Chaos 

3.2 

Ares 

9.5 

12.7 

19.1 

8. Juventae 
Chasma 

3.8 

Maja 

3.3 

7.1 

10.7 

9. Margaritifer 
Chaos 

2.3 

- 

- 

2.3 

3.45 





Total: 296 


* Sediment-water ratio = 0.67 


Table 4.1. Outflow sources used in MODFLOW models, and their associated 
valleys, with the approximate volume of material removed from each. 
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The regional surface recharge is a free parameter and is applied at a rate of 2 x 10' 
10 m/s (6.3 mm.yr' 1 ) for all models, corresponding to estimated south polar ice cap 
basal melting rates of C&P. Recharge elevation coincides with the Martian surface, a 
consequence of local cryosphere melting due either to the insulation of an overlying 
ice sheet, or the proximity of volcanic heat sources (see Discussion). Excess recharge 
is free to drain from the model. Water is also allowed to exit through those parts of 
the upper model boundary corresponding to the circum-Chryse outflow channel 
sources. We assume that discharge occurs because of a disruption in the cryosphere 
and that water is consequently released from the aquifer at an elevation corresponding 
to the pre-disruption cryosphere base. We model the largest nine outflow channel 
sources: Juventae Chasma, Echus Chasma, Margaritifer Chaos, Iani Chaos, Aram 
Chaos, Hydaspis Chaos, Hydraotes Chaos, Aureum-Arsinoes Chaos, and Aromatum 
Chaos (Table 4.1, Figure 4.1). The volume of material removed from each source 
region and its associated outflow channel(s) is taken from Carr (1996). The volume of 
water required to transport this material is determined by the sediment-volume ratio. 
Komar (1980) applied to the large outflow channels a set of well-known terrestrial 
criteria for determining the grain size of different types of stream load. He concluded 
that the outflow channel morphology, coupled with the effects of the reduced Martian 
gravity on water and grain-settling velocities, produced sediment loads in the region 
of 40% by volume (equivalent to a volumetric sediment-water ratio of 0.67), similar 
to values reported for rivers in the semi-arid southwestern United States (Beverage 
and Culbertson, 1964). We thus adopt a sediment-water ratio of 0.67 (Table 4.1), but 
consider it an upper limit. Later we present scaling relationships that estimate outflow 
discharge given smaller values. 
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Figure 4.4 Model topography cast in an equal-area sinusoidal projection spanning 
longitudes from 180° W to 0°, and all latitudes, (a) South Pole and (b) Tharsis 
recharge regions are shown in black and modeled outflow channel sources are shown 
in white. 


We adopt a conservative approach in determining discharge requirements, 
modeling the main chaos regions but only two Chasmata, specifically Echus and 
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Juventae, from which sediment was carried along single obvious channel pathways 
(Kasei and Maja Valles, respectively). We do not include water required to remove 
material from the other large Chasmata since it is unclear what proportions of this 
water were distributed to particular outflow channels and through the aquifer between 
chasmata. We thus account for about half (3 x 10 6 km’. Table 4.1) the total discharge 
estimated to have passed through the circum-Chryse outflow channels, a small 
decrease considering the orders of magnitude uncertainty in sediment- water ratio. 

All outflow sources are assumed to begin discharging simultaneously at the start of 
the model. This simplification, while not entirely physically realistic, measures the 
ability of the recharge region to provide the necessary discharge to a fully active 
Chiyse outflow region. It also provides some insight into the interaction between 
outflow sources and how water is shared between them. 

4.2.3 Recharge and Thermal Properties 

We examine three models, each with a different combination of thermal properties 
and recharge location. In our nominal South Pole recharge model (model 1) we adopt 
the C&P cryosphere thickness derived from an estimate of 95 mW.m' 2 for the 
Hesperian heat flux (Figure 4.2), which in turn was derived from a parameterized- 
convection global thermal model (Stevenson et al., 1983). Annually averaged surface 
temperatures correspond to insolation at the current obliquity and vary from 154 K at 
the pole to 210 K at the equator. Porosity varies with depth according to Equation 4.4, 
with a surface value of 20% (C&P). Surface recharge occurs everywhere within 750 
km of the south pole (C&P, Figure 4.1) and since basal melting precludes the 
presence of a cryosphere beneath the ice cap, recharge is applied to the crustal 
surface. 



87 


Q. 

Q 


to 

m 


o 

£> 

CJ 



-90 -60 -30 0 30 60 90 

S N 


Figure 4.5 Latitudinal dependence of cryosphere base depth for the hydrologic and 
thermal conditions adopted in models 1, 2 and 3. The region of zero cryosphere 
thickness in models 1 and 3 represents basal melting of a polar ice cap, which 
requires the melting isotherm to reach the surface. In model 2, zero cryosphere 
thickness (whose lateral extent in this plot corresponds to a longitude of 245°) is 
achieved by basal melting or volcanic heating (see Discussion). The effects of the 
adiabatic lapse rate in models 2 and 3 are not shown here, since they depend on 
topography. The lower model boundary corresponds to a depth of 8.5 km in all 
models. 


Our Tharsis Recharge model (model 2) tests the ability of Tharsis recharge to 
provide outflow channel discharge, including more conservative assumptions than in 
model 1 . Ft covers the same hemisphere as the nominal model, with the same lateral 
discretization. The assumed Hesperian heat flow of 30 mW.m' 2 is the mean value 
derived by McGovern et al. (2002) based on gravity-topography relationships. 
Cryosphere thickness is thus greater than in model 1 (Figure 4.2). Annual average 
surface temperatures are assumed to be the mean of values at 25° and 60° obliquity 
and are thus somewhat lower at the equator and higher at the pole than at present, 
leading to smaller latitudinal variations in cryosphere thickness. We also include an 
adiabatic lapse rate of 2.5 K.krrf 1 (Zurek et al., 1992) making higher elevations 
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slightly cooler. Porosity varies with depth according to Equation 4.5. Recharge is 
applied to Tharsis elevations above 6000 m (Figure 4.1), since the areal coverage of 
the topography above this elevation is approximately the same as the South Pole 
recharge area (about 1.5 x 10 6 km 2 ), allowing a more balanced comparison of model 
results. 

In our alternative South Pole model (model 3), we test South Pole recharge given 
the surface temperatures and geothermal heat flux of the Tharsis Recharge model. 

4.3. Results 

4.3.1 Nominal South Pole Recharge Model 
Instantaneous outflow channel discharges from model 1 (Figure 4.3) show that all 
outflow sources are able to expel the required volume of water by about 403 Myr, the 
lower bound on the duration of outflow activity (Figure 4.4a). Fulfillment of 
discharge requirements appears to support south pole recharge, but comparison with 
an identical model in which no recharge occurs shows almost the same results (gray 
curves in Figure 4.3). In fact, South Pole recharge contributes negligibly to outflow 
discharge: it is principally the high-lying Tharsis province that empties its initial 
groundwater inventory into the outflow channels, and the relatively thin cryosphere 
allows it to do so in the short time observed. Only the high assumed sediment-water 
ratio makes this possible: smaller values would increase the discharge requirement, 
whereupon contributions from South Pole recharge would be essential. However, the 
result that South Pole recharge contributes negligibly to outflow discharge over 400 
Myr at the modeled permeability, holds regardless of sediment-water ratio. 
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Figure 4.6 Instantaneous discharges from the nine outflow channel sources of 
model 1 . Source numbering corresponds to Table 4. 1 . The results of an identical 
model in which no recharge is applied, are shown by the gray curves. 

We note that the individual outflows sources interact strongly: when one source 
completes its required flow and is shut off, its groundwater supply is diverted to the 
remaining active outflows, which experience surges in discharge. Only the Echus 
Chasma source remains active after 50 Myr, and the hydraulic head at 400 Myr 
(Figure 4.4a) reflects the dominance of this source on the flow pattern. 

It is useful to consider the results of this model scaled to different values of some 
key parameters. First, there is great uncertainty in the sediment- water ratio estimates 
that are so fundamental to water-volume calculations. We suggest that this ratio is 
bounded from below by crustal permeability and outflow discharge duration. 
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(a) 



Figure 4.4 Hydraulic head distribution at 400 Myr in (a) model 1 and (b) model 2. 
Outflow activity ceases at approximately this time in both models. Recharge zone and 
outflow sources are shown in black and white, respectively. Scale bars are in km. 

Specifically, we know that the circum-Chryse outflows occurred from the late 
Noachian to the late Hesperian, with an upper bound of about 900 Myr. A preliminary 
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estimate of the lower bound on sediment-water ratio is therefore that value which 
allows the outflows to fulfill their discharge requirements in exactly 900 Myr. The 
analysis can be broadened by repeating the process for a range of viable 
permeabilities. Manning and Ingebritsen (1999) suggest limits of 10' 17 and 10' 14 m 2 
for the mean continental permeability between 1 and 10 km depth, and we scale these 
to the Martian environment (2 x 10' 16 and 2 x 10' 13 m 2 ). By comparison, Harrison and 
Grimm (2002) derived a minimum permeability for the Martian crust of 10‘ l7 -10' 16 m 2 
by comparing the requirements for hydrothermal alteration to numerical models of 
hydrothermal circulation. Rate of discharge can be scaled with complete accuracy 
only by the time-intensive process of rerunning the model for each permeability 
value. We assume instead that model 1 instantaneous discharge rates summed over all 
nine outflows and averaged over the duration of outflow activity (400 Myr in the 
model) apply to every model of the same permeability, regardless of its duration. 
(This approach can be interpreted as an assumption of steady state conditions under 
which instantaneous discharge rates are constant with time and therefore the same 
regardless of model duration.) We assume further that this time-averaged discharge 
rate then scales linearly with permeability in the same way as instantaneous rates, as 
required by Darcy’s Law. 

In Figure 4.5a, the results are shown in a plot of permeability vs. sediment-water 
ratio. The dot represents model 1 - its permeability value (1.0 x 10" 14 m 2 ) as 
calculated by averaging over all aquifer cells in the model corresponds to an outflow 
activity duration of 400 Myr. Only those points within the shaded region are 
physically permissible since they fall below the upper bound on sediment-water ratio 
(0.67), between the bounds on permeability, and they correspond to total discharge 
that is complete within 900 Myr. The minimum sediment-water ratio, corresponding 
to a permeability of 2 x 10' 13 m 2 and a time of 900 Myr, is 0.015, is about two orders 
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of magnitude higher than ratios inferred for the relatively small outflow channels 
Maumee and Vedra Valles (de Hon et al., 2003). 



Figure 4.5 (a) Scaling of model 1 average permeability and sediment-water ratio. 
Bounds discussed in the text limit points to the gray region, (b) Scaling of model 
1 average permeability and outflow activity duration, (c) Same as (a) but for model 2. 
(d) Same as (b) but for model 2. 


Next, we consider the influence of recharge rate. When model 1 is rerun with 
applied recharge rates as low as 10' 12 m.s' 1 , the outflow discharge is the same as that 
of the nominal model with its applied value of 2 x 10' 10 m.s" 1 . Clearly only a small 
fraction of the applied recharge is infiltrating into the aquifer. Although C&P use 2 x 
10' 10 m.s' 1 as a suggested infiltration rate, we see that the aquifer is not capable of 
conducting this much water. The actual infiltration rate will be similar in magnitude 
to the average discharge rate, which in model 1 is about 5 x 10' lj m.s' 1 (calculated 
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from the average flux of 0.23 mV 1 , assuming all outflow drains are active). Model 1 
is thus operating in a regime in which large perturbations in the recharge rate have no 
effect on infiltration - only changes in permeability will change the infiltration rate. 
We thus look at a range of average permeabilities, estimating the duration of outflow 
activity using the same assumption concerning time-averaged discharge rate (and 
therefore infiltration rate) made with our sediment-water ratio calculations. The 
result (Figure 4.6b) is a linear decrease in outflow duration with infiltration rate and 
permeability (the dot again represents model 1). Outflow activity with a duration of 
exactly 900 Myr is achieved with an average permeability of 5 x 10' 15 m 2 (compare 
with the model 1 average value of 1.0 x 10" 14 m 2 ) and an infiltration rate of 2 x 10' 13 
m.s" 1 . A minimum outflow activity duration, 13 Myr, is imposed by the upper bound 
on permeability (2 x 10" lj m 2 ), and corresponds to a recharge rate of 1 x 10" 11 m.s" 1 , 
still less than the applied recharge rate. 



Time [Myr] 


Figure 7.6 Instantaneous discharges from the nine outflow channel sources of 
model 2. Source numbering corresponds to Table 4. 1 . The results of an identical 
model in which no recharge is applied, are shown by the gray curves. 
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4.3.2 Tharsis Recharge Model 

Instantaneous outflow channel discharges from model 2 (Figure 4.6) show that 
Tharsis recharge is capable of producing the required output within 396 Myr. This is 
remarkable given that the average aquifer permeability is about 4 times smaller than 
that of model 1, principally because of the greater model 1 cryosphere thickness. This 
disadvantage is overcome by the proximity and elevation of Tharsis recharge which 
produces high hydraulic head gradients that start early and are sustained indefinitely. 
Nonetheless, some outflows finish later than in the nominal model. Comparison with 
a model in which no recharge is applied (gray curves in Figure 4.6) reveals 
substantially smaller discharges, with one outflow source (Aureum-Arsinoes, Table 
4.1) still active after 900 Myr. The need for recharge is also observed in the total 
amount of water initially available in the aquifer, namely 5.4 x 10 5 km 3 , less than the 
total required outflow channel discharge (Table 4.1). 

The interaction between sources is similar to that observed in model 1, with the 
addition of a gradual increase in all discharges. This occurs as recharge slowly raises 
the hydraulic head of the region (as predicted by C&P for South Pole recharge; 
compare with gray curves, Figure 4.6), increasing the relative drawdown of the 
outflow channels. For model 1 South Pole recharge, this effect is not noticeable 
during the 400 Myr of outflow activity. 

We look again at the lower bound on sediment-water ratio. The minimum 
permissible sediment-water ratio in this model is 0.0038 (Figure 4.5c), corresponding 
to a permeability of 2 x 10‘ ,J m 2 and a time of 900 Myr. This is approximately one 
order of magnitude smaller than the lower bound for model 1 and only one order of 
magnitude higher than ratios for Maumee and Vedra Valles (de Hon et al., 2003). The 
results of Figures 5a and 5c encapsulate quantitatively the advantage of Tharsis 
recharge over South Pole recharge. The former is capable of producing approximately 
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the same time-averaged outflow discharge rate as the latter, even though average 
aquifer permeability is about four times smaller (ratio of spatially averaged 
permeabilities = 1.0 x 10' l4 /2.6x 10* 15 =3.8). 

Considering a range of infiltration rates (Figure 4.5d, note that the ordinate range is 
different to Figure 4.3b) reveals that outflow activity lasting exactly 900 Myr is 
achieved with an average permeability of 1 x 10' 15 m 2 and an infiltration rate of 2 x 
10' b m.s' 1 . The minimum outflow activity duration is 3.2 Myr, corresponding to an 
infiltration rate of 4 x 10*" m.s" 1 , still less than the applied recharge rate. The results 
of Figures 5b and 5d again quantify the advantage of Tharsis recharge, although in 
slightly different terms. With equal infdtration rates for Tharsis and South Pole 
recharge, the duration of outflow activity is approximately the same despite the large 
difference in average permeability. 

Finally, we look at our choice of recharge (precipitation) area. For a variety of 
average infiltration rates and permeabilities, we plot in Figure 4.7a the recharge area 
required to assure the same flux as model 2 (0.24 m 3 .s*'). The lower and upper bounds 
on permeability (2 x 10* 16 m 2 and 2 x 10* 13 m 2 ) place respective upper and lower 
bounds on recharge area, namely 2 x 10 13 m 2 and 2 x 10 10 m 2 . The upper bound 
corresponds to an area almost half that of the tropical band within 30° of the equator 
(vertical dashed line, Figure 4.7a) and therefore somewhat larger than the Tharsis 
rise. We discuss below the conditions under which recharge might infiltrate the 
aquifer only through some fraction of the applied area. 

4.3.3 Alternative South Pole Recharge Model 

In the first 500 Myr of model 3, the outflow sources do not (with the exception of 
the large Aureum-Arsinoes chaos region. Table 4.1) expel much more water than in 
an identical model with no recharge (Figure 4.8). No source achieves its total required 
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output before 500 Myr, and only two finish before 900 Myr. Again, it is mainly the 
drainage of Tharsis that provides water to the outflows, but with model 2 hydrologic 
and thermal properties this reservoir is not sufficient, as already pointed out in the 
results for model 2. 

For the Aureum-Arsinoes chaos region, the reduced Tharsis inventory means a 
noticeable relative discharge contribution from South Pole recharge (compare the 
black and gray curves for outflow source 3 in Figure 4.8). However, absolute 
discharge is over an order of magnitude less in model 3 than in model 1, and Aureum- 
Arsinoes reaches its discharge requirements in 800, rather than 22 Myr. 

We make some important observations regarding cryosphere thickness. Consider 
model 1 hydrologic properties with an improved Hesperian heat flux estimate of 30 
mW.m'" (McGovern et al., 2002). The resulting south pole cryosphere thickness is 
between 6.5 to 8 km (for melting temperatures ranging from 252 K to 273 K), 
bringing into question the degree of aquifer connectivity. If porosities are low 
enough, it is possible for parts of the cryosphere to occupy the entire uncompacted 
depth of the crust, thereby blocking the lateral flow of groundwater. We calculate 
latitude-dependent aquifer thickness under this scenario for a range of surface 
porosities, in order to quantify aquifer connectivity. 
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Area of Applied Recharge [m 2 ] 



Figure 4.7 (a) Average permeability and infiltration rates required to maintain the 
model 2 flux through a range of applied recharge areas. The dot represents model 2. 
(b) Infiltration rate required to maintain a model 2 outflow activity duration of 396 
Myr given a range of reductions in infiltration area with respect to precipitation area. 
The applied recharge rate of model 2 is indicated by the dotted line. Where the two 
lines cross, all applied recharge infiltrates the aquifer (i.e. no surface runoff). 
Infiltration rates above the dotted line are thus impossible to deliver unless the 
applied recharge is increased. 
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Figure 4.8 Instantaneous discharges from the nine outflow channel sources of model 
3. Source numbering corresponds to Table 4.1. The results of an identical model in 
which no recharge is applied are shown by the gray curves. 


For all surface porosities, the southern hemisphere aquifer is thinnest near the edge 
of the polar cap at about -77° latitude (Figure 4.9; note that no cryosphere is present 
under the south pole ice cap due to basal melting). In the low heat flow model (Figure 
4.9b), surface porosities lower than about 10% result in zero aquifer thickness at this 
latitude, obstructing the northward flow of basal meltwater. This result demonstrates 
the importance of high surface porosity and geothermal heat flux in model 1 (Figure 
4.9a). Hypothetically, if the cryosphere were to grow due to a decrease in surface 
temperatures, the flow of South Pole recharge would be cut off from the outflow 
channels while any equatorial recharge would remain connected. As a caveat, we note 
that the true Martian compaction depth is not known to any great accuracy. If 
compaction occurs at porosities lower than 1%, then South Pole recharge would be 
cut off at higher surface porosities. 
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Figure 4.9 Aquifer thickness as determined by the porosity of Equation 4.4 given a 
variety of surface values n (h and South Pole recharge- In (a) the assumed geothermal 
heat flux is 95 mW.m' 2 (the bold line corresponds to the n 0 of C&P), in (b) it is 30 
mW.m' 2 . In both panels, the aquifer thickness determined by the porosity of Equation 
4.5, with Tharsis recharge, is shown (dotted) for comparison. The indicated lateral 
extent of cryosphere melting for Tharsis recharge corresponds to a longitude of 
245° E. 
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4.4 Discussion 

We elaborate on some problematic aspects of the C&P South Pole recharge model 
and then describe the alternative of Tharsis recharge, which may occur through the 
melting of an ice sheet analogous to South Polar basal melting (C&P), or by direct 
infiltration of precipitation into parts of the aquifer whose cryosphere has melted due 
to volcanic heating. 

4. 4. 1 Problematic aspects of the Nominal South Pole Recharge Model 

The C&P South Pole recharge model requires a groundwater aquifer with a high 
degree of global connectivity, with pressures having to communicate across distances 
as great as 5000 km from south pole to equator. C&P invoke a decrease in outflow 
source elevation with time as evidence of a falling global water table. This argument 
is not supported by the formation process of the Amazonian channels, which may 
have been different from that of the Noachian and Hesperian channels (Car, 1996), 
nor is it supported statistically - the Amazonian channels are small in number and are 
spatially concentrated on the north and northwest flanks of Elysium Mons. With no 
strong evidence of global connectivity, the possibility of regional heterogeneity in 
Martian crustal permeability and porosity remains. Such conditions would drastically 
reduce the already marginal role of South Pole recharge in outflow channel activity 
(as predicted by model 1 ). 

Second, the C&P model has difficulty explaining the provision of water to the 
highest outflow channels, which sometimes exceed the elevation of the south pole 
(Carr, 2002; Coleman et al., 2003). C&P argue that hydraulic head can rise above the 
level of the south pole topography when the polar cap confines the aquifer, but this is 
not physically possible, as the head cannot rise above the maximum recharge 


101 


elevation. C&P's comparison to a recharge well is inappropriate because the latter 
provides external mechanical energy to force water into the ground. The Chryse 
outflows require more water than is initially available in the high elevation Tharsis 
aquifer (given the improved hydrologic and thermal conditions of model 2) and the 
South Pole region is too low to recharge these regions. 

4.4.2 Tharsis Recharge: Ice Sheet with Basal Melting 

As precipitation accumulates to form a growing ice sheet, the melting isotherm that 
defines the lower boundary of the cryosphere rises, releasing melted ice into the 
aquifer. Once the ice sheet reaches a thickness comparable to the pre-existing 
cryosphere, the melting isotherm reaches the ice sheet itself, basal melting begins, 
and proceeds at a rate approaching the precipitation rate. Our Hesperian model of 
Tharsis recharge has an equatorial exosphere thickness, and thus an implied ice sheet 
thickness, of 2 km. If the precipitation rate is equal to the applied recharge rate of 2 x 
10' 10 m.s’ 1 , such a sheet would accumulate in 0.3 Myr, shorter than an obliquity cycle. 
The general climate model of Mischna et al. (2003), which starts from a steady state 
solution at the current obliquity, takes 30 yr to accumulate a total mass of 3 x 10 16 kg 
within 30° of the equator after obliquity is changed to 60°. This accumulation rate (7 
x 10' 10 m.s' 1 ) is linear and, if sufficient water were available (as it would be from 
cyclical outflow discharges during the Hesperian), would build a 2 km ice sheet in 
only 0.09 Myr, roughly the duration of an obliquity cycle. 

Even if deposition rates are not this high, it is likely that ice could remain at the 
tropics during short periods of unfavorable obliquity, allowing accumulation to 
continue during subsequent favorable periods. Mustard et al. (2003) describe 
obliquity cycles with amplitudes as high as 35° that produced a 1.7 Myr ice age that 
ended only 0.4 Myr ago. They argue that ice would continue to accumulate at mid- 
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latitudes throughout this ice age, with only small losses during the traverses through 
low obliquities between cycle maxima. 

Additionally, the presence of a thick (> 3 km) polar ice sheet today implies that 
(non-ice age) sheet growth is not destroyed every 0. 1 Myr, so ice has been preserved 
through unfavorable obliquities by the presence of lag deposits or some other 
mechanism. Indeed, observations of ice-related processes on the flanks of several of 
the Tharsis volcanoes (Tanaka, 1985; Head and Marchant, 2003; McGovern et al., 
2003) and elsewhere at low latitudes (Boynton et al., 2002; Mustard et al., 2001 and 
2003; Christensen, 2003; Greeley and Fagents, 2001) testify to ice preservation 
processes at low latitudes also. A comparison of polar and low-latitude atmospheric 
conditions reveals that tropical summer temperatures at low obliquities would be a 
few 10’s of degrees lower than south polar summer temperatures at high obliquities 
(Ward, 1992; Mischna et al., 2003; Haberle et al., 2003). Retention of ice at the 
tropics during high obliquity is thus expected to be more efficient than retention at the 
poles during low obliquity (Herkenhoff and Murray, 1990; Vasavada et al., 2000). 

4.4.3 Tharsis Recharge: Local cryosphere melting 

There exists an alternative recharge pathway at Tharsis, if the ice sheet is not 
regionally thick enough to induce basal melting. Volcanic activity during the 
Hesperian may have heated significant fractions of the Tharsis crust, resulting in a 
thin cryosphere (Tanaka et al., 2002) that at some points may have disappeared 
altogether (Wilson and Head, 2002). High heat fluxes from local intrusions would 
thus obviate the need for an ice sheet, allowing precipitation to infiltrate (at least on a 
local scale) directly into exposed parts of the aquifer. The disadvantage of this form 
of recharge is that infiltration will not occur at all places experiencing precipitation - 
only some fraction of the precipitation area will experience recharge. In order to 
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places bounds on the scale of this reduction, consider a variant of model 2 in which 
precipitation falls over the same region as before, but infiltrates only through that 
fraction corresponding to the parts with melted cryosphere. The outflows will have to 
remain active for a longer period in order to fulfill their discharge requirements. As 
calculated above, the average infiltration rate that leads to an outflow activity 
duration of 900 Myr is 2 x 10' b m.s' 1 . The ratio of this rate and the model 2 rate (5 x 
1 O' 1 3 m.s' 1 ) thus approximates the lower bound on the ratio of infiltration and 
precipitation areas, namely 0.4. Put differently, infiltration must occur over at least 
40% of the model 2 precipitation area for discharge requirements to be met. 

If instead we demand that the total flux of water into the aquifer remains at a 
constant value of 0.24 nr’.s' 1 (and therefore that outflow activity duration remains at 
396 Myr), we can calculate the precipitation rate required for a given ratio of 
infiltration-precipitation area (Figure 4.7b). The dotted line represents the applied 
recharge of 2 x 10' 10 m.s' 1 , and for points on the solid line above this threshold, 
infiltration cannot be delivered at the indicated rate unless the applied recharge is 
increased. However, by assuming that permeability scales with infiltration as before, 
its upper bound (2 x 10' 13 m 2 ) limits the infiltration rate to 4x 10' 11 m.s' 1 (Figure 
4.5d), which does not exceed the applied recharge of model 2. This maximum 
infiltration rate corresponds to an infiltration-precipitation area ratio of 1.3%, or a 
total area of melted cryosphere only 19600 km 2 . 

4.4.4 Other Outflow Channels 

Finally, we make some remarks concerning outflow channels not belonging to the 
circum-Chryse group. These include Dao, Reull, and Harmakhis Valles on the 
northeast rim of Hellas Planitia; Ma’adim and Mangala Valles in the Amazonis 
Planitia; and the Amazonian channels on the northern and northwestern flanks of 
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Elysium Mons. Only Ma’adim and Mangala Valles do not lie close to volcanic 
centers. Mangala Valles issues from a graben that likely disrupted the local 
cryosphere. Recent work suggests that Ma’adim Vallis formed due to catastrophic 
overflow of a paleo-lake (Irwin et al., 2002) and therefore might not conform to the 
model of groundwater discharge from the cryosphere. 

4.5 Conclusions 

C&P explain the delivery of water to the outflow channels by the communication 
of confined pore fluid pressures across a globally connected aquifer from a South 
Polar recharge zone northwards to the tropics, where water discharges through a 
disruption in the cryosphere and erodes the circum-Chryse outflows. We performed 
numerical simulations to test both this hypothesis and an alternative in which the 
aquifer is recharged over the Tharsis rise, and therefore is not required to be globally 
connected. Our conclusions follow. 

1. In a nominal model with hydrologic and thermal properties assumed by C&P, 
South Pole recharge contributes little to outflow discharge. The assumed permeability 
results in fulfillment of discharge requirements only because the high sediment-water 
ratio allows the entire discharge volume to be supplied by water drained from the 
Tharsis aquifer. At lower sediment-water ratios, the weak contribution of South Pole 
recharge will not be able to produce required discharges in the allotted time unless 
permeabilities are increased. 

2. In a model with improved, yet more conservative properties, Tharsis recharge is 
able to meet discharge requirements in approximately the same time as the nominal 
South Pole model. In this case, recharge is an essential contributor to outflow 
discharge, asserting a strong, sustained influence. 
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3. South Pole recharge with the improved properties of the Tharsis recharge model 
performs worse than the nominal model, with only 2 of the 9 modeled outflows 
completing their discharge requirements in 900 Myr. 

Simple linear scaling of our Tharsis recharge model results indicates that for 
average permeabilities as high as 2 x 10' b m 2 , the minimum sediment-water ratio that 
allows discharge requirements to be met in the maximum allowed time of 900 Myr, is 

0.0038. Alternatively, with a high sediment-water ratio (0.67) and high permeability, 
the region over which recharge occurs could be as small as 2 x 10 4 km 2 . If the 
recharge region remains at the assumed value of 1.5 x 10 6 km 2 , high permeabilities 
might allow discharge requirements to be met in only 3.2 Myr. 

We argue that recharge through Tharsis occurs at high obliquity when ice is 
perennially stable in the tropics. Infiltration of water into the aquifer may take place 
in at least two ways: 

1 . Basal melting of an ice sheet. In analogy to the melting of a South Polar ice cap 
(C&P), we suggest that the insulating effects of an accumulating ice sheet might raise 
the melting isotherm to intersect with the base of the sheet, resulting in the steady 
infiltration of meltwater into the Tharsis aquifer, from whence it flows down a steep 
hydraulic head gradient to the outflow channels. 

2. Local cryosphere melting. The presence of volcanic activity in Tharsis 
throughout the Hesperian suggests that the melting isotherm may have locally and 
temporarily risen to the surface, allowing the direct infiltration of precipitation. 
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CHAPTER 5 
CONCLUSIONS 

we conclude this work with a summary of essential results, a brief analysis of key 
unsolved issues influencing my conclusions, and a description of factors necessary to 
improve the accuracy of my calculations. 

5.1 Valley Networks 

Hydrotherm discharge attending crustal formation processes may have been 
sufficient to supply the erosional requirements of the ancient Southern Highlands 
valley networks that correlate spatially with crustal magnetic anomalies. This result 
is supported by numerical calculations of groundwater convection due to the intrusion 
of a magma chamber. The nominal model of a 50 km 3 chamber in host rock of 
permeability 10‘ 16 m 2 provides enough discharge, when building that part of the crust 
covered by the valley networks, to fall within the bounds of the water volume 
required for valley erosion. Larger chambers, possibly more appropriate to the 
martian environment, produce even more discharge than the nominal model. Factors 
producing a reduction in discharge include evaporation and, more importantly, water 
table elevation. 

Discharge requirements place a lower bound on host rock permeability of 10‘ 17 
m 2 , while alteration requirements suggest a minimum of 10' 16 m 2 . We note that 
discharge requirements depend on the sediment-water ratio, which for the valley 
networks might range from 10‘ 3 to 0.25. 
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5.2 Long Runout Landslides 

Numerical simulations of martian landslides with a variety of material rheologies 
suggest that long runout landslides in the Valles Marineris may have resulted from 
the fluidization of the falling debris due to pore fluid pressure. A fluidization 
rheology which partitions failure into an initial frictional stage and a final power law 
stage, produces landslide runouts that match observed profiles in five of seven 
martian landslides modeled. Frictional sliding applied to the entire duration of failure 
only reproduces long runouts with unrealisitically high pore pressures and, along with 
other rheologies such as Bingham flow, cannot reproduce certain aspects of the final 
profile, such as the toe height and rear displacement. 

Pore pressures required by our fluidization rheology are high, but physically 
tenable. The equatorial location of the landslides implies that cryosphere thickness 
may have been low enough to allow the stable presence of pore water in the failure 
zone, which in the largest landslides is a few kilometres below the surface. 

Simulations of an Olympus Mons aureole indicate that it is unlikely to be the 
product of catastrophic landsliding since unrealistically high pore pressures are 
required for all rheologies. Formation processes suggested by other workers do, 
however, require the action of water or ice. 

5.3 Outflow Channels 

Clifford and Parker (200 1 ) suggest that the confined Hesperian aquifer was 
recharged at the South Pole in order to stimulate global groundwater leading to 
outflow channel formation. My suggestion that recharge over Tharsis at high 
obliquity is a more efficient means of driving outflow discharge is strongly supported 
by numerical simulations that compare the two recharge scenarios. For the chosen 
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crustal permeability, South Pole recharge has a weak effect on outflow discharge 
during the supposed period of formation (at most 900 Myr), whereas Tharsis recharge 
dominates discharge production throughout the duration of outflow activity. 

As with valley networks, the discharge requirements depend on sediment-water 
ratio, which for our models is assumed to be at a maximum value of 0.67. Smaller 
sediment-water ratios require a greater volume of discharge through the outflow 
channels, which can be accommodated in the allowed time if permeability is 
increased. Simple scaling relations show that a reasonable maximum (average) 
permeability of 2 x 10' 13 m 2 allows a sediment-water ratio as low as 0.4% for outflow 
discharge fed by Tharsis recharge. 


5.4 Historical Relationship Between Landforms and Environment 

The relative timing of valley networks, long runout landslides, and outflow 
channels is shown in Figure 5.1, along with crustal heat fluxes and cryosphere 
thicknesses assumed by C&P and myself (Chapter 4). 

There is still some debate regarding Noachian atmospheric conditions. Surface 
runoff implied by valley network formation seems to suggest temperatures above 
freezing, but younger Amazonian valleys, such as those on Alba Patera, probably 
formed under cold atmospheric conditions similar to those of the present. If the 
Noachian atmosphere was indeed cold, an explanation of the sapping morphology of 
the valley networks must incorporate the effects of a cryosphere on the order of 1 km 
thick (Figure 5.1). As shown in Chapter 2, hydrothermal discharge is capable of 
melting through such a cryosphere with a significant, but not crippling reduction in 
surface discharge. Warm hydrothermal discharge might remain in the liquid phase 
long enough to cause significant erosion, as suggested by the Alba Patera valleys. The 
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Figure 5.1 (a) Relative timing of valley networks, outflow channels, and long runout 
landslides. Dashed lines indicate decreased production of the indicated landform. (b) 
Crustal heat fluxes assumed by C&P, and calculated by McGovern et al. (2002). The 
latter curve represents approximately the lowest values calculated by McGovern et al. 
(2002). Hesperian values correspond to those used in models 1 and 2 of chapter 4. (c) 
Equatorial cryosphere thicknesses result from the crustal heat fluxes of (b), and from 
surface temperatures corresponding to an obliquity of 25° (“C&P”), and to an average 
of obliquities between 25° and 60° (“Tharsis Recharge Model”). The linear increase 
over the N-H epoch is a response to the depletion of an assumed thick, warm 
Noachian atmosphere. A thin, cold Noachian atmosphere gives rise to the dashed 
curves. 


argument for hydrothermal ly driven discharge is in fact stronger under cold 
atmospheric conditions, since alternatives involving precipitation, which would be in 
solid form, are unlikely to produce runoff. At best, precipitation would contribute to 
valley network formation by infiltrating into locally melted parts of the cryosphere (in 
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analogy to Tharsis recharge during the Hesperian) and mixing with hydrothermal ly 
circulating groundwater. 

Figure 5.1 shows a decline in valley network production during the late Noachian 
and early Hesperian. Ostensibly this decline is due to the development of a 
cryosphere following the depletion of a warm Noachian atmosphere (Carr, 1996). If 
the Noachian atmosphere was cold and dry, a decline in valley network formation 
might still occur in response to waning crustal formation (and therefore hydrothermal 
activity), although it may not be as rapid. 

Uncertainty in the role of groundwater arises from apparently incompatible 
inferences drawn from mineral assemblages in the Valles Marineris. An analysis of 
infrared MGS TES data by Staid et al. (2001) indicates the presence of unaltered 
olivine, implying persistently dry conditions. The deposit is not mentioned explicitly 
by other workers (e.g. Hoefen and Clark, 2001; Clark et al., 2002) but other deposits, 
mostly in ancient terrain (e.g. Terra Meridiani, Syrtis Major, southern Acidalia 
Planitia, and Terra Tyrrhena; Hoefen and Clark, 2001) give rise to the same 
conclusion: a dry early Mars (Christensen et al., 2001a). Also detected in the Valles 
Marineris are deposits of crystalline gray hematite, a product of aqueous alteration. 
The hematite is closely associated with the interior layered deposits (Christensen et 
al., 2001b) and implies the long-term presence of liquid water on a local scale. 
Groundwater flow in the Valles Marineris is also indicated by morphological aspects 
of the interior deposits, and by long runout landslides, and its influence is suggested 
indirectly by inter-chasma groundwater features (Carr, 1996) and the proximity of the 
outflow channels. 

An apparently paradoxical picture thus emerges of continuously dry conditions 
coexisting with regions of local water stability. Further analysis of the olivine 
deposits, especially in the Valles Marineris, is required, and must take into account 



the degree of alteration, which is difficult to detect with the available infrared data 
(Kraft et al., 2002). 

With regard to long runout landslides, note that the equatorial cryosphere 
thickness (Figure 5.1) is never greater than 2 km, supporting my suggestion that the 
failure- zones of Valles Marineris landslides may have been saturated with water. 


5.5 Future Improvements 

The acceleration in our understanding of water on Mars, as alluded to in the 
Introduction, is a continually developing process. There are still several gaps in our 
knowledge of Martian aquifer conditions that in this work occasionally translate into 
uncertainties in quantitative results. Groundwater simulations cannot accurately 
model a real physical environment if little is known about the nature of its 
permeability and porosity. To date, hydrologic properties of the Martian aquifer have 
been based on extrapolations from terrestrial data, which are themselves incomplete 
and difficult to interpret. Further measurements of the Martian environment would 
increase our understanding of both the porosity structure (e.g. seismic data) and the 
present location of ice and groundwater (e.g. radar). 

The weak constraints on sediment-water ratio in valley network and outflow 

channel formation also need to be addressed. Detailed analysis of channel and 
bedform morphologies, coupled with improved knowledge of the size, distribution, 
and history of the water inventory of Mars, would tighten these constraints. 
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